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ABSTRACT 

A  coupled  ice-ocean  numerical  model  is  developed  which 
improves  the  simulation  of  the  annual  cycle  and  interannual 
variations  in  ice  cover  in  the  Arctic.  The  model  is  a  further 
development  of  the  work  by  Semtner  (1987) .  Although  the 
accuracy  of  the  simulated  ice  concentration  is  increased,  the 
annual  cycle  of  ice  coverage  is  still  exaggerated.  Several 
experiments  are  conducted  to  determine  the  importance  of 
incorporating  a  fully  interactive  ocean,  to  select  an  optimum 
strength  parameter  for  use  in  the  ice  rheology,  to  investigate 
the  model's  sensitivity  to  changes  in  the  albedo  of  the  frozen 
surface  and  to  determine  the  relative  importance  of  the 
various  dynamic  and  thermodynamic  forcing  mechanisms.  The 
regional  dependence  of  these  mechanisms  and  an  assessment  of 
two  statistical  analysis  techniques  used  to  measure  model 
improvement  are  also  examined. 

Inclusion  of  a  fully  prognostic  ocean  component  vice  a 
ten-year  mean  ocean  cycle  in  the  model  improves  the 
correlation  of  simulated  ice  concentration  fields  with 
observed  data.  This  is  the  case  for  all  regions  in  the 
Arctic;  for  both  the  annual  cycle  and  interannual  variations 
of  the  ice  cover.  A  reduced  strength  parameter  value, 
P  =hxlO^,  is  found  to  improve  the  simulation  of  the  ice 
thickness  distribution  with  increased  overall  thickness  and 
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better  compression  north  of  the  Canadian  Archipelago  and 
Greenland.  In  contrast  to  results  using  ice  models  without 
a  fully  prognostic  ocean  component,  this  model  is  quite 
insensitive  to  changes  in  the  frozen  surface  albedo. 
Exceptions  are  evident  where  the  ocean  heat  flux  into  the 
mixed  layer  is  small  and  the  ice  is  thin. 

At  the  spatial  (110  km)  and  temporal  (monthly)  scales  used 
here,  the  heat  provided  by  the  ocean  appears  to  be  the 
dominant  mechanism  controlling  the  position  of  the  ice  edge 
and  the  extent  of  the  ice  pack.  Within  the  pack,  it  is  the 
dynamic  forcing  and,  in  particular,  the  wind  forcing  which 
controls  the  ice  thickness  and  thickness  distribution.  The 
ocean  circulation  below  the  mixed  layer  appears  to  position 
the  heat  underneath  the  MIZ.  The  MIZ  is  also  the  region  where 
the  ice  thickness  tends  to  decrease  through  divergence.  The 
linkage  between  the  subsurface  heat  and  the  thinned  ice  cover 
is  apparently  controlled  by  conditions  at  the  surface  and  the 
resulting  response  of  the  mixed  layer. 
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I.   INTRODUCTION 

Interest  in  the  complex  geophysical  processes  occurring 
in,  above  and  around  the  Arctic  Ocean  continues  to  increase 
as  the  economic,  strategic  and  political  values  of  this  region 
become  more  important.  In  the  case  of  economic  concerns,  vast 
oil  reserves  are  being  exploited  in  this  area  and  the  search 
for  more  reserves  continues.  Much  of  Western  Europe,  northern 
Canada  and  the  USSR  are  supplied  via  shipping  routes  which  are 
dependent  on  weather  and  ice  conditions  in  the  Arctic  regions. 
Arctic  waters  also  contain  some  of  the  most  valuable  fishing 
grounds  in  the  world;  however  our  ability  to  use  them  is  again 
dependent  on  the  state  of  ice,  water  and  weather. 

The  Arctic  Ocean  spans  the  shortest  distance  between  North 
America  and  the  USSR  and  is  therefore  of  obvious  strategic 
importance.  However  the  extreme  environmental  conditions 
encountered  there  make  military  operations  very  difficult  and 
expensive.  Naval  vessels  used  in  northern  waters  must  be 
properly  designed  if  they  are  to  function  with  any 
effectiveness  and  their  performance  will  normally  be 
significantly  degraded.  The  reguirement  for  an  ice  breaking 
capability  is  obvious.  In  addition,  the  highly  variable 
temperatures,  surfaces,  and  material  compositions  of  water, 
ice  and  snow  have  a  dramatic  effect  on  both  the  air  and  water 
boundary  layers.    This  in  turn  affects  electro-magnetic. 


electro-optical,  and  acoustic  transmission  and  reception. 
Sea-ice  is  a  strong  source  of  ambient  acoustic  noise  and  radar 
scatter.  Detection  of  surface  or  low  level  air  targets  by 
radar  or  submarines  by  passive  or  active  acoustic  means 
becomes  very  difficult. 

Several  countries  border  the  Arctic  mediterranean  and  have 
strong  political  motives  for  establishing  or  extending  their 
sovereignty  and/or  economic  control  as  far  offshore  as 
possible.  The  presence  of  several  "chokepoints"  within  the 
region  provides  the  opportunity  to  influence  operations  in  the 
entire  Arctic  Basin  by  controlling  only  limited  areas. 

Obviously,  a  comprehensive  knowledge  of  the  Arctic 
environment,  including  the  ocean,  atmosphere  and  sea-ice,  and 
a  capability  to  predict  this  environment  with  some  accuracy 
would  be  of  great  value.  The  expense  of  establishing  a 
complete  observational  network  is  prohibitive.  An  important 
role  must  be  played  by  numerical  modelling  which  has  the 
potential  to  provide  both  a  full  analysis  of  the  current 
conditions  based  on  limited  data,  and  a  prediction  capability 
for  future  conditions. 

A.   NUMERICAL  MODELLING  PERSPECTIVE 

It  is  recognized  that  the  atmosphere,  ice  and  ocean 
interact  in  a  highly  complex  nonlinear  fashion.  Consideration 
of  one  of  these  three  components  cannot  be  realistically  done 
in  total  isolation  of  the  other  two.  However  a  complete 


understanding  of  all  the  physics  involved  in  describing  each 
parameter  and  its  interaction  with  the  others  is  lacking. 
Therefore  our  ability  to  mathematically  describe  the  physics 
involved  is  also  lacking.  Furthermore,  when  significant 
theoretical  advances  have  been  made  in  the  past,  memory  and 
speed  limitations  in  computing  facilities  have  been  an 
obstacle.  Even  with  today's  supercomputers,  compromises  must 
be  made  to  balance  computing  requirements  (memory,  expense  and 
time)  against  the  accuracy  desired.  In  order  to  minimize 
these  compromises,  the  three  component  problem  has  generally 
been  broken  up  into  its  three  individual  components  of 
atmosphere,  ocean  and  ice.  The  two  components  which  are  not 
being  specifically  examined  in  a  model  are  then  parameterized 
in  some  fashion  or  they  are  represented  by  prescribed 
empirical  data. 

Numerical  modelling  was  first  applied  to  the  atmosphere. 
Impetus  was  provided  by  public  as  well  as  military  interests. 
Great  strides  have  been  made  by  meteorologists  in  predicting 
the  weather  by  numerical  models.  Numerous  atmospheric 
numerical  models  of  varying  complexity  have  been  developed  and 
are  currently  in  use  by  the  various  National  Weather  Centers. 
However,  despite  a  great  deal  of  effort  and  expense, 
deterministic  prediction  of  daily  weather  fluctuations  is 
presently  limited  to  five  to  ten  days.  Theoretical  arguments 
suggest  that  the  ultimate  practical  range  of  such  predictions 
is  about  two  weeks. 


Numerical  modelling  of  the  global  oceans  requires  a 
thousandfold  increase  in  computer  power.  Compared  to 
numerical  modelling  of  the  global  atmosphere,  the  horizontal 
resolution  of  the  global  oceans  must  be  ten  times  as  great, 
and  the  numerical  time  step  therefore  ten  times  as  small,  to 
accommodate  the  smaller  radius  of  deformations  and  faster 
gravity  wave  speeds.  Ocean  modelling  is  further  complicated 
by  irregular  borders,  islands,  widely  varying  density 
structures,  internal  and  boundary  mixing  processes,  complex 
and  varied  boundary  layers  and  a  marked  paucity  of 
observations.  Nevertheless  with  the  major  contributions  of 
Bryan  (1969)  and  numerous  others  such  as  Takano  (1974), 
Semtner  (1974),  Cox  (1984)  and  Semtner  and  Chervin  (1988), 
viable  three-dimensional,  baroclinic  numerical  models  have 
been  developed  which  simulate  the  state  of  the  world's  oceans 
quite  well.  These  models  are  currently  in  wide  use  in 
circulation  studies,  tracer  distribution  analysis,  etc. 

Numerical  modelling  of  the  ice  has  also  progressed 
significantly  over  the  past  twenty  years.  The  physics  of  this 
problem  are  now  considered  to  be  fairly  well  understood. 
Prediction  of  the  ice  both  in  areal  extent  and  thickness 
distribution,  is  known  to  be  a  function  of  dynamic  forcing 
(both  internal  and  external)  and  thermodynamic  forcing.  The 
work  of  Hibler  (1979,  1980)  has  provided  details  of  the 
important  dynamic  processes.  Maykut  and  Untersteiner  (1969, 
1971)    have    constructed    a    detailed    one-dimensional 


thermodynamic  model  which  provided  a  comprehensive  view  of  the 
thermodynamic  aspects.  In  the  short  term  (days) ,  wind  stress 
appears  to  be  the  dominant  factor  in  ice  prediction.  On  this 
time  scale  the  thermodynamic  effects  of  freezing  and  melting 
and  the  dynamic  effects  of  weak  ocean  currents  are  small.  As 
a  result,  predictions  of  sea-ice  on  this  time  scale  can  be 
reasonably  accurate  with  only  minimal  consideration  (simple 
parameterizations)  of  the  ocean's  influence.  However, 
attempts  to  predict  ice  on  a  seasonal  basis  using  either  a 
thermodynamic  or  a  dynamic  model,  or  some  combination  of  the 
two,  without  adequate  consideration  of  the  ocean's  influence, 
have  met  with  only  limited  success.  The  ice  and  ocean  are 
closely  linked  over  the  long  time  scale  and  the  interactions 
which  occur  between  them  must  be  accounted  for. 

The  U.S.  Navy's  current  operational  ice  forecasting  model 
run  at  the  Fleet  Numerical  Oceanography  Center  (FNOC)  is 
called  the  NORDA/FNOC  Polar  Ice  Prediction  System  (PIPS) .  One 
of  the  limitations  of  this  model  is  its  use  of  the  monthly 
average  ocean  currents  and  ocean  heat  fluxes  produced  from  the 
Hibler  and  Bryan  (1987)  ice-ocean  model.  This  model  was 
forced  for  several  years  of  integration  by  repeatedly  using 
a  single  year  of  observed  atmospheric  forcing  data  (December 
1978  to  November  1979) .  This  year  was  chosen  because  it  was 
the  "FGGE"  year  in  which  a  large  number  of  drifting  buoys  were 
in  the  Arctic  Basin  which  could  provide  a  check  for  the  model 
currents.   The  time  period  was  not  chosen  for  its  similarity 


to  an  average  condition  and,  as  noted  by  Hibler  and  Bryan 
(1987) ,  it  was  in  fact  substantially  different  from  the  long 
term  average.  Additionally,  the  ocean  below  the  mixed  layer 
was  constrained,  with  a  three  year  damping  time,  to  the 
observed  circulation.  Therefore,  the  ocean  fluxes  and 
currents  calculated  from  the  output  of  this  model  may  very 
well  be  limiting  the  accuracy  of  the  PIPS  model  predictions. 
This  problem  has  yet  to  be  resolved. 

The  present  work  is  aimed  at  improving  our  understanding 
and  capability  to  numerically  predict  the  location  and 
thickness  of  Arctic  ice  on  a  seasonal  time  scale.  This  will 
be  done  through  the  further  development  of  a  linked  ice-ocean 
numerical  model  using  observed  atmospheric  forcing  and  a  fully 
prognostic  ocean  model  formulation.  Several  sensitivity 
studies  will  also  be  conducted  to  gain  a  better  appreciation 
of  the  relative  importance  of  the  various  mechanisms 
controlling  the  ice  cover  and  to  improve  their  representation 
in  the  model.  The  remainder  of  this  chapter  briefly  describes 
previous  efforts  in  ice-ocean  modelling  which  have  provided 
the  background  for  this  work.  The  specific  goals  and 
objectives  of  this  research  are  also  listed.  Chapter  II 
provides  the  geophysical  context  in  which  the  model  operates 
and  Chapter  III  is  a  description  of  the  model  itself  including 
the  assumptions,  boundary  conditions  and  forcing  fields  under 
which  it  is  run.  Chapter  IV  examines  the  importance  of 
including  a  fully  interactive  ocean  model  in  studies  of 


interannual  ice  cover  variations  and  Chapter  V  examines  the 
effect  of  ice  strength  on  those  same  variations  and  the  ice 
thickness  distribution.  Chapter  VI  investigates  the 
sensitivity  of  the  model  to  frozen  surface  albedo  changes  and 
Chapter  VII  is  a  comprehensive  examination  of  the  various 
dynamic  and  thermodynamic  processes  included  in  the  model  and 
their  relative  importance  to  the  evolution  of  the  ice  cover. 
Chapter  VIII  addresses  the  value  of  comparing  the  model 
simulations  statistically  and  Chapter  IX  is  a  final  summary 
which  includes  the  conclusions  from  this  work. 

B.   PREVIOUS  ARCTIC  MODELLING 

Coupled  ice-ocean  models  of  the  Arctic  on  seasonal  time 
scales  were  developed  only  recently.  Construction  of  such 
models  was  initiated  primarily  in  response  to  the  lack  of 
quantified  success  of  previous  large  scale  ice  modelling 
efforts  by  numerous  investigators.  General  seasonal  trends 
and  basin-wide  distributions  similar  to  observed  fields  were 
achieved  in  these  earlier  efforts.  However,  the  accuracy  of 
regional  ice  distributions  and/or  the  phase  of  the  annual 
cycle  were  often  poor. 

Manabe  and  Stouffer  (1980)  developed  a  coupled  ice-ocean- 
atmosphere  model  which  used  a  one-dimensional  mixed  layer 
ocean  and  a  motionless  slab  with  no  leads  for  the  ice.  The 
atmospheric  model  provided  daily  variable  forcing.  This 
produced  a  wintertime  ice  cover  about  double  that  normally 


observed  in  the  Arctic.  Washington  et  al.  (1976,  1980)  used 
a  simple  thermodynamic  ice  model  but  did  not  include  any  ice 
dynamics.  Their  model  produced  an  excess  of  ice  in  the  North 
Atlantic  in  all  but  the  summer  months  and  ice  thicknesses  in 
the  central  Arctic  were  less  than  half  those  interpreted  from 
submarine  surveys  (LeSchack  et  al.,  1971;  LeSchack,  1980). 

Once  the  emphasis  shifted  to  prediction  of  the  ice,  the 
atmospheric  portion  of  the  coupled  models  was  generally 
dropped  in  favor  of  prescribed  atmospheric  forcing  and  more 
elaborate  ocean  and  ice  models.  The  detailed  thermodynamic 
ice  model  of  Maykut  and  Untersteiner  (1969,  1971)  included  the 
effects  of  ice  salinity,  brine  pockets,  shortwave  radiation 
heating,  vertical  density  variations,  conductivity  and 
specific  heat  of  ice  and  water.  However,  this  model  required 
a  great  deal  of  computation  to  reach  an  equilibrium  state. 
This  made  it  impractical  for  application  to  a  large  scale 
three-dimensional  coupled  model.  Semtner  (1976a)  simplified 
this  complex  model  by  parameterizing  several  of  the  processes. 
With  the  use  of  a  three-layer  version  he  was  able  to  closely 
reproduce  the  results  of  Maykut  and  Untersteiner  (1969,1971) 
in  a  wide  variety  of  simulations,  yet  this  simplified  model 
had  fewer  computational  requirements  than  one  layer  of  ocean 
in  a  multi-level  ocean  model.  This  made  it  ideal  for 
inclusion  into  a  large  three-dimensional  gridded  coupled 
model . 


The  Parkinson  and  Washington  (1979)  ice  model  used  a 
further  simplified  "zero-layer"  thermodynamic  ice  model  (also 
proposed  and  tested  successfully  by  Semtner  (1976a)),  in 
conjunction  with  a  dynamic  forcing  model.  The  dynamics  were 
based  on  stresses  from  wind,  water,  Coriolis  force,  internal 
ice  resistance  and  ocean  surface  tilt.  The  atmosphere  was 
represented  by  monthly  climatological  fields  of  temperature 
and  pressure  and  the  ocean  was  represented  by  a  simple  3  0  m 
deep  ocean  mixed  layer  with  a  constant  vertical  heat  flux. 
This  model  produced  a  reasonable  yearly  cycle  of  sea  ice 
extent;  however,  the  ice  thickness  distribution  did  not 
compare  well  with  observational  evidence.  The  combined 
thermodynamic-dynamic  ice  model  concept  was  further  improved 
by  Hibler  (1979,  1980).  He  also  used  the  "zero-layer" 
thermodynamic  approach  of  Semtner  (1976a)  but  combined  it  with 
a  more  complex  dynamic  model  and  used  daily  observed 
atmospheric  forcing  data.  The  ocean  was  again  represented  by 
a  simple  fixed-depth  mixed  layer  with  a  constant  vertical  heat 
flux  applied  beneath  it.  The  essential  difference  in  the 
Hibler  approach  was  to  link  the  dynamics  of  ice  motion  to  the 
ice  thickness  by  allowing  the  ice  interaction  to  become 
stronger  as  the  ice  grew  in  thickness  and/or  concentration. 
In  order  to  do  this  consistently,  a  viscous-plastic  ice 
rheology  was  used.  Arctic  simulations  with  this  model  showed 
that  the  spatial  distribution  of  ice  thickness  was  strongly 
influenced  by  the   ice  dynamics.     The   results   compared 


favorably  with  the  ice  distribution  observations  of  LeSchack 
(1980)  and  were  a  significant  improvement  over  previous 
efforts.  The  improved  ice  dynamics  as  well  as  daily 
atmospheric  forcing  permitted  a  much  better  representation  of 
lead  formation.  The  proper  representation  of  leads  was 
considered  important  because  heat  and  salt  fluxes,  as  well  as 
ice  formation  rates  are  nominally  an  order  of  magnitude  larger 
in  open  leads  than  under  the  pack  ice. 

Further  analysis  of  the  Semtner  (1976a)  thermodynamic  ice 
model  was  conducted  by  Semtner  (1984a) .  He  concluded  that  an 
oversimplification  of  the  thermodynamics  tends  to  exaggerate 
the  seasonal  sea-ice  cycle.  The  "three-layer"  version  of  the 
thermodynamic  model  involves  twice  as  much  computation  as  the 
"zero-layer"  version  but  does  provide  a  much  better 
representation  of  the  ice.  The  computational  requirements 
remain  negligible  compared  to  that  of  the  ocean  model. 
Therefore  the  "three-layer"  model  was  recommended  for  use  in 
large  coupled  models  where  accurate  representation  of  the 
regional  and  temporal  changes  in  the  ice  field  are  important. 

Several  experiments  have  been  conducted  to  determine  the 
importance  of  atmospheric  forcing.  Three  years  of 
interannually  variable  forcing  was  applied  by  Hibler  and  Walsh 
(1982)  to  the  Hibler  (1979)  ice  model,  but  omitting  inter- 
active ocean  circulation  effects.  Their  results  again 
demonstrated  the  importance  of  ice  dynamics  for  achieving 
accurate  ice  forecasts.  Quantitative  agreement  with  observed 
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interannually  varying  ice  volume  and  transport  was  achieved; 
however,  a  noticeable  seasonal  bias  was  evident.  The  ice 
extent  and  thickness  in  winter  and  spring  for  all  years  was 
excessive  while  the  extent  and  thickness  of  summer  ice  was  too 
small  with  too  much  melting  and  excessive  open  water.  They 
concluded  that  accurate  parameterization  of  oceanic  effects 
and  a  better  treatment  of  the  snow  cover  would  probably  be 
necessary  in  order  for  Arctic  ice  to  be  accurately  modeled. 
Walsh  et  al.  (1984,  1985)  increased  the  interannually  varying 
atmospheric  forcing  period  to  30  years  and  used  a  more 
elaborate  thermodynamic  model  with  the  Hibler  (1979)  ice 
dynamics.  The  seasonal  bias  in  ice  concentration  was  again 
evident  in  this  model  and  it  also  showed  considerable 
sensitivity  to  the  number  of  vertical  layers  in  the 
thermodynamic  portion. 

An  Arctic  ocean  circulation  model  was  developed  by 
Semtner  (1976b)  as  an  extension  to  the  numerical  ocean  model 
of  Bryan  and  Cox  (Bryan,  1969)  .  Many  of  the  observed  features 
of  the  Arctic  Ocean  and  Greenland  Sea  circulation  were 
reproduced  by  this  Arctic  model  despite  the  simplifications 
of  mean  annual  wind  forcing  and  a  non-interactive  ice  cover. 
In  order  to  develop  an  ice  model  suitable  for  seasonal  ice 
predictions  and  daily  variable  atmospheric  forcing,  Hibler  and 
Bryan  (1984,  1986)  coupled  the  Semtner  (1976b)  Arctic  ocean 
model  with  a  somewhat  simplified  Hibler  (1979)  thermodynamic- 
dynamic  ice  model.    This  resulted  in  the  first  seasonal, 
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coupled,  ice-ocean  numerical  model.  A  linear  term  was 
included  in  the  ocean  model  which  damped  the  ocean's 
temperature  and  salinity  values  below  the  mixed  layer  to 
climatology.  This  "robust-diagnostic"  forcing  was  chosen  to 
have  a  three  year  relaxation  time.  Multi-year  drift  from 
climatology  was  prevented;  yet  shorter  term  variations, 
particularly  in  the  mixed  layer,  could  still  be  prognostically 
simulated.  The  inclusion  of  the  oceanic  component  greatly 
improved  the  simulated  ice  edge,  removing  much  of  the  seasonal 
biases  evident  in  the  earlier  work  by  Hibler  and  Walsh  (1982) 
and  Walsh  et  al.  (1984,  1985).  Hibler  and  Bryan  (1987)  showed 
that  many  of  the  features  of  ice  growth  and  decay,  in  various 
regions  of  the  Arctic,  result  from  a  combination  of  numerous 
influences,  often  with  opposing  effects.  For  example,  the 
East  Greenland  ice  edge  is  forced  southward  by  the  East 
Greenland  Current  and  by  winter  freezing.  However,  due  to  a 
strong  vertical  heat  flux  from  the  ocean,  the  amount  of  ice 
melted  actually  exceeds  the  amount  frozen  in  this  region.  The 
southward  extension  of  the  ice  zone  is  thus  significantly 
inhibited. 

The  Hibler  and  Bryan  (1984,  1986)  model  simulated  the 
seasonal  variations  in  ice  cover  more  accurately  than  any  of 
the  previous  modelling  efforts.  However,  Semtner  (1987)  noted 
that  the  three  year  diagnostic  constraint  may  have  dominated 
the  forcing  for  the  ocean  model  and  thereby  controlled  the  ice 
simulation.   The  constraint  also  prevented  ocean  circulation 
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and  ice  phenomena  occurring  on  multi-year  time  scales  (bottom 
water  formation,  carbon  dioxide  effects,  etc.)  from  being 
forecast.  Additionally,  the  climatological  ocean  was 
unrealistic  in  the  variable-current,  shallow-water  regions  of 
the  Barents  and  East  Siberian  Seas.  A  fully  prognostic, 
coupled,  ice-ocean  model  would  have  no  such  limitations  and 
was  therefore  developed  (Semtner,  1987).  This  model, 
hereafter  referred  to  as  SM,  is  the  basis  of  the  current  work. 
SM  incorporated  the  best  features  of  the  previous  modelling 
efforts.  The  Semtner  (1976b)  ocean  model,  optimized  for  high 
speed  computers,  was  combined  with  the  "three-layer" 
thermodynamic  ice  model  (Semtner,  1976a)  and  a  modified 
Hibler  (1979)  dynamic  model.  The  daily  forcing  ice  dynamics 
formulation  of  Hibler  (1979)  required  more  computation  than 
a  complete  ocean  model.  This  was  not  considered  an  equitable 
distribution  of  computation  load.  The  Hibler  dynamics  were 
modified  to  use  monthly  averaged  forcing  instead  of  daily 
forcing,  yet  still  retain  the  first-order  dynamic  effects 
(Hibler,  1988) .  This  also  permitted  a  much  longer  time  step 
and  longer  integrations.  A  simple  30  year  mean  seasonal  cycle 
was  computed  from  the  forcing  data  used  in  Walsh  et  al.  (1985) 
for  atmospheric  forcing.  Monthly  values  of  water  inflow  at 
the  appropriate  salinity,  temperature  and  volume  from  rivers 
and  the  southern  ocean  boundaries  were  specified.  These  water 
mass  fluxes  were  invariate  in  time.  The  model  produced 
realistic  representations  of  the  average  seasonal  ice  cycle 
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in  distribution  and  concentration.  The  ocean  circulation, 
including  the  major  surface  currents,  water  masses,  and 
possible  water  transformation  regions  were  also  well  simulated 
to  the  limits  of  the  resolution.  Results  from  SM  emphasized 
the  high  degree  of  thermodynamic  and  dynamic  interaction 
between  the  ice  and  ocean. 

C.   OBJECTIVES 

The  objectives  of  this  work  were  to  gain  further 
understanding  of  the  relative  importance  of  mechanisms 
controlling  the  ice  cover  in  the  Arctic  Basin  and  to  develop 
a  numerical  model  capable  of  accurately  simulating  the 
operationally  and  climatologically  important  interannual 
variations  of  ice  cover  in  this  region.  In  particular,  the 
following  points  unique  to  this  work  were  to  be  investigated: 

-  Determine  the  importance  of  including  a  fully  prognostic 
and  interactive  ocean  component  in  a  linked  ice-ocean 
model.  The  importance  was  to  be  assessed  relative  to  the 
accuracy  of  simulations  of  interannual  variability  in  the 
Arctic  ice  cover. 

-  Determine  the  sensitivity  of  the  model  to  changes  in  the 
ice  rheology  strength  parameter  and  subsequently  select 
an  optimum  value  of  this  parameter  for  incorporation  into 
the  model . 

-  Determine  the  sensitivity  of  the  model  to  changes  of  snow 
ind  ice  albedo  and  if  appropriate  incorporate  a  more 
sophisticated  albedo  representation  into  the  model. 

-  Determine  the  relative  importance  of  the  various  dynamic 
and  thermodynamic  ice  cover  forcing  mechanisms  included 
in  the  model,  particularly  the  two  features  inherent  in 
the  ocean  model  (vertical  ocean  heat  flux  and  ocean 
current  stress) . 
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Assess  the  regional  dependency  of  the  above  ice  cover 
control  mechanisms. 

Assess  the  statistical  significance  of  the  differences 
between  the  various  model  run  results. 


D.   METHODOLOGY 

This  research  was  accomplished  through  the  utilization  and 
further  development  of  the  coupled  numerical  model  of  Semtner 
(1987) .  The  Semtner  (1987)  model  is  the  most  freely 
prognostic,  physically  representative  and  computationally 
efficient  of  the  large-scale,  linked,  ice-ocean  numerical 
models  previously  developed  for  the  Arctic.  Considerable 
skill  has  been  demonstrated  with  the  model  in  reproducing  the 
annual  cycle  of  ice  extent  with  climatological  forcing.  Most 
importantly,  the  inclusion  of  an  interactive  ocean  reduced  the 
pronounced  tendency  of  previous  ice-only  simulations  to 
exaggerate  the  annual  ice  cycle. 

Numerous  model  integrations  were  conducted.  In  each 
experiment  the  model  was  modified  to  examine  some  particular 
mechanism  which  was  expected  to  influence  the  interannual 
variability  of  the  ice  cover.  In  the  majority  of  cases  the 
model  was  integrated  for  ten  years.  However,  some  of  the 
sensitivity  testing  required  a  large  number  of  repeated  runs 
and  in  those  cases  the  model  runs  were  limited  to  one 
integration  year.  Additionally,  in  the  sensitivity  runs  where 
inclusion  of  an  interactive  ocean  was  of  limited  concern,  the 
model  was  run  with  a  prescribed  ten-year  average  ocean 
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condition  vice  a  prognostic  ocean  component.  The  ice  cover 
simulations  were  evaluated  from  January  1971  to  December  1980 
inclusive.  These  years  were  chosen  because  the  observed  ice 
concentration  data  set,  to  which  the  modelled  ice 
concentrations  were  compared,  was  most  accurate  in  this 
period.  This  was  especially  true  after  1972  due  to  the 
assimilation  of  microwave  remotely  sensed  data.  Parkinson  et 
al.  (1987)  was  also  an  excellent  source  of  comparison  data  for 
the  years  1973-1976. 

The  observed  ice  concentration  fields  were  provided  by 
Walsh,  and  were  the  same  as  those  used  in  Fleming  (1987)  but 
converted  through  bi-linear  interpolation  to  the  grid  used 
here.  The  observed  ice  concentration  fields  were  regarded  as 
the  "true"  monthly  ice  concentrations.  However,  they  were  in 
fact  the  average  of  data  received  from  several  sources,  each 
with  its  own  source  of  error  (Walsh  and  Johnson,  1979) . 
Nevertheless,  these  fields  were  still  the  most  comprehensive 
and  accurate  representations  of  real  data  available  and 
covered  a  long  enough  time  period  to  permit  viable  statistical 
analysis.  All  120  months  of  observed  ice  concentration  were 
contoured  to  gain  a  qualitative  appreciation  of  the  annual 
cycle,  its  variability  and  absolute  extent.  This  also 
provided  the  baseline  view  for  visual  comparison  with  the 
simulated  ice  concentration  data. 

The  seasonal  average  ice  draft  contours  from  Bourke  and 
Garrett  (1987)  were  the  primary  source  of  comparison  data  for 
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the  simulated  ice  thicknesses  and  thickness  distributions. 
Although  these  contours  represent  "observed"  data  from 
submarine  transits,  the  data  were  heavily  averaged  in  time  and 
space  to  produce  the  contours.  Nevertheless,  the  general 
thickness  and  distribution  patterns  did  provide  an  indication 
of  the  true  seasonal  ice  thickness  distribution. 

The  monthly  varying  atmospheric  forcing  fields  were  also 
provided  by  Walsh.  These  fields  were  produced  from  the  same 
data  used  by  Walsh  et  al.  (1985)  and  were  also  the  source  of 
the  3  0-year  mean  forcing  using  in  Semtner  (1987) . 

It  was  shown  by  Fleming  (1987)  and  Garcia  (1988)  that 
correlations  between  ice  concentration  and  thermodynamic 
forcing  variables  (sea  surface  temperature  (SST)  and  air 
temperature)  and  dynamic  forcing  variables  (winds  and  previous 
ice  condition)  vary  significantly  between  regions  which 
experience  different  dynamic  ocean  conditions.  To  gain  some 
appreciation  for  these  regional  differences  but  still  maintain 
the  amount  of  data  manipulation  within  reasonable  limits,  the 
complete  grid  space  was  divided  into  four  regions  (see  Figure 
1.1).  The  divisions  were  selected  based  on  the  general  ocean 
current  regime.   The  four  regimes  represented  were: 

-  Region  1     X=l,18      Beaufort  Gyre. 

-  Region  2      X=19,27     Transpolar  Drift,  Ellesmere  Is. 

and  Greenland  ice  convergence 
zone. 
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Figure  1.1 


Geography  of  the  Arctic  Ocean  as  represented  in 
the  model.  The  small  boxes  indicate  the  centers 
of  land  gridpoints.  Latitude  circles  every  two 
degrees  are  shown  over  the  ocean  gridpoints.  The 
arrows  point  to  the  mouths  of  the  indicated  rivers 
and  the  numbers  indicate  the  regions  (from 
Semtner,  1987) . 
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-  Region  3      X=28,34     Transpolar  Drift  (TPD) ,  East 

Greeenland  Current  (EGC) . 

-  Region  4     X=3  5,45     Norwegian  Current. 

where  X  defines  the  location  of  the  region  on  the  X  axis  of 
the  grid.  Each  region  extends  from  1  to  42  on  the  Y  axis  of 
the  grid.  In  those  cases  where  the  modified  model  simulations 
produced  ice  fields  more  representative  of  the  observed  ice 
fields  than  before  the  modifications  were  made,  the 
appropriate  changes  were  incorporated  into  the  model. 
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II.   ARCTIC  CLIMATOLOGY 

This  chapter  presents  the  geophysical  characteristics  of 
the  Arctic  region  emphasizing  those  features  which  are  most 
relevant  to  numerical  modelling  of  the  ice  cover  and  ocean. 

The  primary  dynamic  driving  forces  for  sea  ice  growth  and 
movement  are  the  air  stress  resulting  from  wind-induced 
surface  drag  and  the  water  stress  resulting  from  ocean 
currents.  Variations  in  wind  forcing  are  determined  by  the 
passage  of  synoptic  scale  weather  systems  having  spatial 
scales  of  100  to  1000  km  and  durations  of  several  days.  The 
response  of  the  sea  ice  to  mesoscale  features  of  the 
atmosphere  and  ocean  is  to  produce  eddies  and  other 
fluctuations  with  length  scales  of  approximately  10  km  and 
time  scales  of  hours  to  days.  Inertial  oscillations  of 
drifting  sea  ice  have  similar  periods.  Longer  term 
fluctuations,  with  time  scales  of  weeks  to  seasons,  represent 
responses  to  atmospheric  and  ocean  forcing  integrated  over 
equivalent  time  scales.  No  firm  evidence  exists  to  suggest 
that  small  time  and  spatial  forcing  is  unimportant  to  the 
seasonal  evolution  of  the  ice  cover.  However  numerical 
modelling  of  the  Arctic  ice  has  not  developed  sufficiently  at 
present  to  resolve  such  forcing  for  long  periods  over  the 
entire  Arctic  basin.  Therefore,  the  emphasis  in  this  work 
will  be  on  the  monthly  time  and  larger  spatial  scale  features. 
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When  examining  fluctuations  varying  on  time  scales  of  a 
season  to  several  decades,  the  dominant  signal  in  sea  ice 
variability  is  clearly  the  annual  cycle  (Parkinson  et  al., 
1987).  Figure  2.1  shows  the  25-year  (1953-1977)  envelopes  of 
the  positions  of  the  Arctic  ice  edge  at  the  end  of  August  and 
February.  At  most  longitudes  the  seasonal  change  from  the 
summer  to  the  winter  ice  edge  positions  is  considerably 
greater  than  the  25-year  range  of  extremes  for  a  particular 
month.  The  overall  annual  freeze-melt  cycle  is  primarily  a 
response  to  the  seasonally  varying  insolation,  albeit  with  a 
large  amount  of  feedback  and  interaction  with  the  ocean. 

A.   ATMOSPHERE 

The  climate  in  the  Arctic  region  is  largely  determined  by 
four  aspects  (Sater  et  al.,  1971): 

-  A  distinctive  cycle  of  prolonged  periods  of  daylight  with 
net  positive  surface  heat  balance  followed  by  prolonged 
periods  of  darkness  with  a  net  negative  surface  heat 
balance. 

-  A  very  persistent,  cold-cored,  circumpolar,  upper  level 
vortex  which  steers  the  mesoscale  surface-level  weather 
systems. 

-  A  surface  cover  of  snow  or  ice  for  periods  of  the  year, 
which  causes  significant  changes  in  surface  albedo  and 
therefore  surface  heat  absorption. 

-  A  strong  temperature  inversion  above  snow  or  ice  surfaces 
which  acts  to  decouple  the  surface  from  higher  level 
atmospheric  activity. 

Many  factors  are  involved  in  determining  the  net  surface 

heat  balance  for  any  particular  region  of  the  Arctic.   These 

will  not  all  be  examined  here.   Suffice  it  to  say  that  the 
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Figure  2 . 1 


Composite  maximum  and  minimum  ice  extents  at  the 
end  of  August  and  February.  The  hatched  region 
is  the  difference  between  extreme  positions  of  the 
0.5  ice  concentration  lines  based  on  25  August 
grids  and  25  February  grids  (1953-1977)  .  The  dots 
represent  the  grid  points  used  in  digitization  of 
monthly  ice  concentration  data  (from  Walsh  and 
Johnson,  1979)  . 
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main  factor  is  believed  to  be  the  insolation  rate.  Figure  2.2 
summarizes  the  effects  in  a  radiation  balance  diagram 
determined  for  the  Barrow  Strait  region.  The  diagram  shows 
long  periods  of  net  heat  loss  which  correspond  to  the  polar 
winter  and  periods  of  large  ice  production.  The  long  heat 
loss  periods  are  followed  by  short  intense  periods  of  net  heat 
gain  corresponding  to  the  long  daylight  hours  and  the  ice 
melting  of  summer.  This  figure  is  reasonably  representative 
of  the  heat  balance  in  the  Arctic  in  general.  However,  a  term 
which  is  missing  and  which  is  regionally  and  spatially  quite 
variable  is  the  advective  heat  flux,  Q^.  This  term  represents 
the  heat  supplied  to  the  surface  from  vertical  and  horizontal 
advection  in  the  ocean.  Mixing  processes  could  also  be 
included.  Huyer  and  Barber  (1970)  suggested  that  this  term 
was  negligible  in  the  Barrow  Strait.  However,  in  other 
regions  of  the  Arctic,  large  differences  in  latitudinal  extent 
of  the  ice  exist  as  can  be  seen  in  Figure  2.1.  This  suggests 
that  Qy  may  be  very  significant  to  the  surface  heat  balance  in 
other  regions  in  the  Arctic. 

Surface  weather  systems,  especially  cyclones,  are 
primarily  generated  in  regions  of  strong  differential  heating. 
Significant  differences  in  surface  temperatures  occur  at  land- 
sea,  ice-sea,  and  warm-cold  current  boundaries.  Several  such 
boundaries  are  found  in  the  Arctic  in  fairly  constant 
positions.   Weather  systems  are   therefore  more  prevalent  in 
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Figure  2.2 


The  radiation  balance  showing  the  short-wave 
incident  radiation,  Q,  the  absorbed  short-wave 
radiation,  Q^,  the  effective  back  radiation,  Q^^, 
and  the  net  radiation,  Q^  (from  Huyer  and  Barber, 
1970) . 
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these  regions,  but  their  subsequent  influence  on  other  areas 
is  dependent  on  the  upper  level  steering  and  the  conditions 
of  the  surface  over  which  they  pass.  The  cyclones  will  tend 
to  increase  in  intensity  if  they  remain  over  relatively  warm 
open  water,  which  provides  heat  and  moisture.  However,  they 
will  tend  to  weaken  and  occlude  if  they  travel  over  cold  and 
dry  (ice  or  land)  surfaces.  These  factors  imply  that 
considerable  variability  exists  in  the  dynamic  atmospheric 
forcing  between  different  regions  and  in  different  seasons  in 
the  Arctic.  Climatological  averages  or  even  monthly  averages 
of  atmospheric  forcing  will  undoubtedly  mask  some  of  this 
variability. 

Variations  in  surface  albedo  cause  the  percentage  of  solar 
radiation  absorbed  to  differ  by  as  much  as  70%  between 
different  types  of  snow,  ice  and  open  water.  Robinson  et  al., 
(1986)  have  used  comprehensive  sets  of  satellite-derived  data 
and  ground  confirmation  data  to  identify  four  ice  surface 
albedo  classes  in  the  Arctic.  Albedos  range  from  0.80  for  a 
surface  with  fresh  snow  covering  95%  of  the  ice  to  0.29  for 
a  heavily  ponded  surface  with  less  that  10%  bare  snow  or  ice. 
Kukla  and  Robinson  (1980)  determined  that  0.10  was  a 
reasonable  albedo  for  open  water  in  the  high  latitudes.  These 
differences  cause  the  melting  and  freezing  processes  to 
accelerate  rapidly  once  started.  Changes  in  the  ice  cover  can 
be  very  dramatic  over  short  periods  of  time.  The  Arctic 
stratus  deck,  common  to  the  region  in  summer,  also  influences 
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surface  albedo  significantly,  but  this  effect  is  very 
difficult  to  quantify.  Considerable  variability  exists  in  the 
literature  as  to  the  correct  albedo  to  use  in  large  scale 
modeling.  Numerous  authors  have  shown  that  numerical 
thermodynamic  ice  models  can  be  very  sensitive  to  this 
variable  (e.g.,  Shine  and  Henderson-Sellers,  1985;  Ross  and 
Walsh,  1987;  Semtner,  1976a).  Inclusion  of  a  realistic  albedo 
representation  is  therefore  difficult  but  probably  important 
in  an  accurate  numerical  sea-ice  model. 

Open  leads  probably  account  for  less  than  10%  of  the 
surface  area  encompassed  by  the  Arctic  ice;  however  the  heat 
and  moisture  exchanges  which  occur  above  them  can  be  an  order 
of  magnitude  higher  due  to  the  large  differences  in 
temperature  and  radiation  characteristics  between  seawater  and 
the  pack  ice  surface.  This  results  in  localized  areas  of  high 
static  instability,  cloud  formation  and  large  heat  and 
moisture  fluxes.  Leads  are  generally  too  small  to  resolve 
using  current  large-scale  numerical  ice  models;  however,  they 
are  important  enough  to  require  that  their  effects  be 
parameterized . 

The  strong  temperature  inversion  common  to  much  of  the 
Arctic  is  the  result  of  the  very  cold  ice  and  land  surfaces. 
In  winter  the  inversion  layer  can  deepen  to  the  850  mb  level, 
becoming  most  intense  under  calm,  clear  anti-cyclonic 
conditions  (Sater  et  al.,  1971).  The  inversion  is  a  very 
persistent  feature  which  will  quickly  reform  after  periods  of 
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strong  winds,  cloud  cover  or  storm  precipitation  as  long  as 
surface  temperatures  remain  relatively  low.  In  summer  this 
feature  is  confined  to  the  polar  ice  cap  as  surface 
temperatures  over  the  ice-free  water  at  more  southern 
latitudes  become  too  high.  The  inversion  acts  to  insulate  or 
decouple  the  surface  from  upper  level  winds  and  also  tends  to 
put  a  cap  on  any  moisture  which  may  be  absorbed  in  the 
boundary  layer.  The  result  is  relatively  calm  (average  less 
than  5  m/s)  surface  wind  conditions  and  a  predominance  of 
stratus  clouds,  especially  when  the  melt  season  commences. 
Care  must  be  taken  to  ensure  that  surface  winds  calculated 
from  higher  level  pressure  charts,  if  used  as  forcing  fields, 
take  this  feature  into  account  to  ensure  realistic  wind 
forcing  on  the  ice. 

B.   WATER  MASSES 

The  waters  of  the  Arctic  Seas  are  often  described  on  the 
basis  of  temperature  and  salinity.  As  such,  they  are 
comprised  of  three  main  water  masses:  Arctic  Surface  Water, 
Intermediate  or  Atlantic  Water,  and  Deep  or  Bottom  Water 
(Coachman  and  Aagaard,  1974).  Arctic  surface  water  is 
generally  limited  in  depth  to  about  200  m.  It  has  the  most 
variable  characteristics  and  can  be  modified  by  the  weather, 
the  season,  and/or  the  physical  environment.  Temperatures  in 
this  layer  vary  from  -1°C  to  over  2°C.  The  salinity  may  be 
uniform  to  approximately  50  m,  below  which  a  sharp  halocline 
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increases  the  salinity  to  about  34.5  ppt  at  the  bottom  of  the 

layer.   The  variety  of  conditions  in  the  surface  water  is 

evident  in  Figure  2.3,  where  vertical  profiles  of  temperature 

and  salinity  from  a  variety  of  Arctic  basins  are  plotted. 

Coachman  and  Aagaard  (1974,  p.  9)  state  that: 

The   most   important   processes   conditioning   and 
modifying  the  surface  layer  are: 

-  Addition   of   mass   (fresh  water)   from   the   land, 
primarily  from  the  large  Siberian  rivers; 

-  Addition  of  fresh  water  locally  through  melting  of 
ice; 

-  Heat  gain  through  absorption  of  solar  radiation  in 
non-ice-covered  areas  during  summer; 

-  Concentration  of  salt  and  hence  increase  of  density 
of  surface  water,  through  freezing  of  ice; 

-  Heat  loss  to  the  atmosphere  through  any  open  water 
surface,  including  leads  in  the  central  Arctic  pack 
ice;  and 

-  Inflow  and  subsequent  mixing  of  Atlantic  and  Pacific 
waters. 

Processes  1,  2,  and  3  normally  occur  only  from  June  to 

September  and  lead  to  a  decrease  in  water  density.   These 

buoyant  waters  form  a  surface  cap  which  absorbs  radiation  and 

warms.   Therefore,  in  summer,  ice  free  regions  tend  to  have 

warmer  and  less  saline  surface  layers.  In  areas  where  the  ice 

does  not  recede,  surface  temperatures  remain  near  freezing  as 

incoming  energy  is  used  to  melt  the  ice,  but  surface  layer 

salinities  are  reduced  due  to  ice  melting. 
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Figure  2 . 3 


Vertical  profiles  of  temperature  and  salinity  for 
various  northern  high-latitude  basins  (from 
Coachman  and  Aagaard,  1974) . 
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Processes  4  and  5  have  the  greatest  impact  in  winter.  In 
some  areas,  such  as  the  shelf  waters,  conditions  are  such  that 
the  water  becomes  dense  enough  to  penetrate  into  the 
intermediate  Atlantic  layer  or  in  extreme  cases  even  form 
bottom  water  (Aagaard  et  al.,  1985).  However,  in  general,  the 
strong  pycnocline  at  the  base  of  the  surface  layer  prevents 
mixing  due  to  surface  density  changes  from  penetrating  below 
200  m  (Coachman  and  Aagaard,  1974) . 

As  a  consequence  of  the  halocline-derived  mixing  barrier, 
process  6  does  not  exert  a  great  influence  on  the  surface 
water  variability.  However,  regions  of  high  surface 
salinities  (33  to  34.5  ppt)  are  found  in  the  Greenland  and 
Labrador  Seas  and  east  Baffin  Bay.  These  high  salinities  are 
due  to  the  advection  of  North  Atlantic  surface  water  into  the 
Arctic  via  warm  surface  currents  from  the  south.  One  further 
consequence  of  the  halocline  barrier  is  that  the  ice  cover  is 
often  insulated  from  the  relatively  warm  Atlantic  layer. 
However,  in  some  locations  like  the  continental  slopes,  the 
Atlantic  layer  is  forced  to  shallower-than-normal  depths. 
Vigorous  surface  mixing  can  then  break  through  the  halocline 
and  vertical  heat  fluxes  can  occur  (Coachman  and  Aagaard, 
1974)  .  This  results  in  the  ice  melting  from  the  bottom  or  not 
forming  at  all,  increasing  the  area  of  open  water. 

Many  mixing  events  are  small  in  scale  but  complicated  and 
intense.  In  a  large  scale  ice  model  these  effects  should  be 
included,  but  as  was  the  case  for  leads,  they  need  to  be 
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parameterized  in  some  manner  due  to  the  scale  limitations  of 
the  model . 

The  second  water  mass  is  called  the  Atlantic  layer.  This 
warm  and  salty  (35.0  to  35.1  ppt)  water  originates  as  North 
Atlantic  water  flowing  into  the  Arctic  basin  through  the  Fram 
Strait  and  the  Barents  Sea  (Weigel,  1987).  It  extends  from 
200  m  to  900  m  with  temperatures  above  0°C.  A  temperature 
maximum  of  approximately  0.4  5°C,  observed  throughout  the 
central  and  western  Arctic  Basin,  occurs  between  3  00  m  and  500 
m  depth,  dependent  on  location.  The  salinity  gradually 
decreases  to  approximately  34.9  ppt  in  the  Arctic  Ocean  and 
Greenland  Sea  and  approximately  34.6  ppt  in  Baffin  Bay.  The 
intermediate  layer  is  a  source  of  heat  which  can  be  tapped  by 
vertical  circulation.  As  a  result  it  may  play  an  important 
role  in  the  net  heat  balance  of  the  surface  mixed  layer.  The 
net  heat  balance  in  the  mixed  layer  ultimately  determines 
whether  or  not  ice  will  form.  The  ocean  model  used  here  has 
been  shown  to  simulate  this  warm  intermediate  layer  quite  well 
(Semtner,  1984,  see  Figure  2.4). 

The  third  water  mass,  the  deep  and  bottom  waters,  have 
temperatures  below  0°C.  The  salinity  is  nearly  constant  from 
the  bottom  of  the  Atlantic  layer  to  the  ocean  floor.  The 
intermediate  Atlantic  water  and  deep  water  are  advected  into 
and  out  of  the  Arctic  seas  from  adjacent  areas,  principally 
through  Fram  Strait  in  the  North  Atlantic.   Both  water  masses 
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Figure  2.4  Vertical  sections  of  observed  temperature  and 
salinity  along  the  prime  meridian,  from  Coachman 
and  Aagaard  (1974)  (upper  panels) ;  and  vertical 
sections  of  simulated  temperature  and  salinity 
across  the  model  domain  at  approximately  20 
degrees  E  (lower  panels) .  The  lower  panels  have 
a  stretched  vertical  scale,  (from  Semtner,  1984b) 
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are  nearly  isohaline  and  therefore  of  nearly  uniform  potential 
density. 

A  temperature  difference  exists  between  the  deep  and 
bottom  waters  of  the  Canadian  and  Eurasian  Basins.  The  deep 
Canadian  Basin  averages  approximately  -0.45°C  while  the 
Eurasian  basin  temperature  is  approximately  -0.79°C  (Aagaard 
et  al.,  1985).  The  deep  and  bottom  waters  of  the  two  basins 
are  kept  separated  by  the  Lomonosov  Ridge  which  acts  as  a  sill 
between  them. 

Observations  of  the  vertical  structure  at  various  stations 
have  been  taken  over  long  time  periods  and  in  many  different 
locations.  A  remarkable  similarity  in  the  profiles  has  led 
to  the  conclusion  that  the  Arctic  basins  are  in  a  long-term 
dynamic  steady-state  condition  (Coachman  and  Barnes,  1961) . 
It  was  further  noted  that  observed  distributions  of  Arctic 
water  properties  are  a  result  of  continuing  processes  within 
the  basins.  Therefore,  surface  water  temperature-salinity  (T- 
S)  profiles  reflect  the  local  modifying  processes,  while  T-S 
profiles  for  depths  below  200  m  reflect  the  common  origin  of 
the  water. 

The  net  heat  balance  ultimately  determines  how  much  ice 
is  produced  or  melted.  The  heat  flux  provided  by  the  ocean 
is  an  important  parameter  in  determining  the  net  heat  balance 
at  each  grid  point  as  noted  in  Chapter  I.  The  ocean  component 
of  a  linked  ice-ocean  model  should  therefore  have  sufficient 
resolution  to  reasonably  define  the  boundaries  of  the  various 
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water  masses  in  three  dimensions.  It  must  also  be  able  to 
simulate  the  water  motion  which  advects  oceanic  heat  to  the 
ice  and  convects  the  cold  saline  water  to  the  bottom. 

C.   CURRENTS 

The  surface  circulation  of  the  Arctic  Basin  has  been 
derived  from  satellite  observations,  the  trajectory  of  ice 
islands,  buoys  and  floe  stations,  and  the  analysis  of  ship's 
tracks.  The  circulation  of  the  Arctic  waters  is  due  both  to 
water  density  differences  and  wind  forcing.  Large  anomalies 
in  the  flow  (compared  to  the  long-term  mean  currents)  are 
often  observed;  however,  the  long-term  mean  surface  currents 
have  been  reasonably  well  documented  in  the  northern  Atlantic 
and  Arctic  Seas  (Krauss,  1986).  Figure  2.5  covers  the  area 
of  interest  in  this  study  and  indicates  the  major  ocean 
currents.  Recent  moored  current  measurements  have  allowed  a 
reassessment  of  the  boundary  and  interior  circulation  in  the 
Arctic  Basin  (Aagaard,  1988) .  Aagaard  emphasizes  the 
important  role  of  mesoscale  eddies  in  the  Canadian  Basin  and 
proposes  a  separate  cyclonic  flow  for  the  Eurasian  Basin. 
This  newest  assessment  is  very  similar  to  the  simulated 
Atlantic  layer  circulation  in  Semtner  (1987)  and  is  a 
significant  change  from  the  previous  basin-wide  cyclonic 
circulation  proposed  by  Coachman  and  Aagaard  (1974) .  Figure 
2.6  shows  this  comparison. 
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Figure  2.5  Schematic  of  the  large-scale  horizontal 
circulation  patterns  in  the  surface  waters  of  the 
Arctic  region  (from  Parkinson  et  al.,  1987). 
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Figure    2 . 6 


Comparison  of  the  simulated  circulation  in  the 
Atlantic   layer   (565   m)   (upper  diagram) (from 
Semtner,   1987)  with  the  reassessed  subsurface 
circulation  presented  by  Aagaard  (lower  diagram) 
(from  Aagaard,  1988) . 


36 


The  Arctic  Ocean  contains  two  main  surface  flow  patterns. 
In  the  Canadian  Basin,  the  Beaufort  Sea  gyre  rotates 
clockwise.  It  is  driven  primarily  by  the  mean  wind.  Across 
the  Lomonosov  Ridge,  the  Transpolar  Drift  Stream  (TPD)  flows 
directly  across  the  Eurasian  Basin,  from  the  Laptev  Shelf  to 
the  Fram  Strait  where  it  exits  as  the  East  Greenland  Current 
(EGC) .  Several  small,  weak  gyres  split  off  the  Transpolar 
Drift  Stream  as  it  passes  close  to  the  various  islands  off  the 
coast  of  the  USSR. 

Circulation  and  ice  movement  patterns  in  the  peripheral 
Atlantic  seas  of  the  Arctic  Ocean  are  dominated  by  two  major 
currents  systems,  one  warm,  the  other  cold  (Krauss,  1986). 
The  cold  currents  support  ice  growth  and  enhance  the  spatial 
extent  of  sea  ice,  while  the  warm  currents  advect  heat  and 
melt  the  ice  or  preclude  its  formation.  The  southerly  extent 
of  the  ice  edge  between  regions  of  dissimilar  temperatures  can 
vary  as  much  as  30  degrees  of  latitude. 

Sea  ice  off  the  east  coast  of  Greenland  is  found  to 
originate  primarily  in  the  Arctic  basin  and  as  a  result  can 
be  quite  thick.  The  EGC  continually  carries  a  wide  belt  of 
Arctic  pack  ice  southward  through  Fram  Strait.  As  the  ice 
season  develops,  the  ice  edge  extends  further  south  down  the 
east  Greenland  coast.  In  extreme  years  the  ice  edge  will 
extend  as  far  south  as  Kap  Farvel  or  as  far  north  as  70 
degrees  N  (Sater  et  al.,  1971). 
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strong  contrasts  in  ice  coverage  are  also  observed  within 
the  Baffin  Bay-Davis  Strait  region.  The  warm  West  Greenland 
Current  (WGC)  follows  the  bathymetric  contours  of  the  southern 
and  western  Greenland  continental  shelf,  keeping  the  southwest 
coast  of  Greenland  ice  free  during  most  winters.  However,  the 
cold  Labrador  Current  contributes  to  the  large  concentrations 
of  ice  found  along  Canada's  east  coast,  well  south  of  Baffin 
Island.  This  current  also  carries  the  majority  of  icebergs 
southward  through  the  Grand  Banks  fishing  zones,  the  Hibernia 
oil  fields  and  into  the  North  Atlantic  sea  lanes  (CIA,  1978) . 

A  reasonably  accurate  representation  of  the  ocean  currents 
is  obviously  important  if  sea  ice  is  to  be  accurately 
simulated.  The  currents  provide  a  direct  influence  through 
dynamic  forcing  as  well  as  a  thermodynamic  influence  via  heat 
advection. 

D.   ICE 

Sea  ice  is  a  feature  which  contributes  to  the  uniqueness 

of  the  oceanography  in  the  Arctic  region.    Coachman  and 

Aagaard  (1974)  state: 

The  general  oceanographic  consequences  of  a  perennial 
or  seasonal  ice  cover  are: 

-  The  water  temperature  of  the  near-surface  layer  in  the 
presence  of  ice  is  always  maintained  close  to  the 
freezing  point  for  its  salinity  by  the  change  of  phase 
process. 

-  Salt  is  excluded  from  the  ice  to  a  varying  extent, 
but  the  water  under  the  ice  is  always  enriched  in  salt 
by  any  ice  growth.  The  dependence  of  water  density 
on  temperature  and  salinity  is  such  that  close  to  the 
freezing  point,  density  is  almost  solely  a  function 
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of  salinity.  Therefore,  ice  formation  can  increase 
the  density  locally  and  some  vertical  convection  may 
result. 

-  In  the  transfer  of  momentum  from  the  atmosphere  to  the 
ocean,  the  wind  must  act  on  the  sea  through  the 
intermediary  of  the  ice. 

The  Arctic  Ocean  ice  pack  is  confined  by  a  nearly 
continuous  boundary  of  land.  The  associated  constraint  on 
equatorward  transport  is  a  major  reason  why  Arctic  ice 
survives  longer  and  develops  into  more  complex  forms  than  ice 
found  in  southern  polar  regions.  An  annual  net  heat  loss  and 
stratification  of  the  underlying  water  also  contribute  to  the 
longevity  of  Arctic  ice. 

Ice  formation  in  open  water  starts  in  autumn.  As  days 
grow  shorter  and  nights  longer,  the  amount  of  solar  insolation 
decreases.  Since  long-wave  radiation  from  the  Arctic  remains 
approximately  constant  during  the  year,  the  energy  budget 
changes  from  a  net  gain  in  summer  through  equilibrium  to  a  net 
loss  in  the  fall  (Figure  2.2).  The  relatively  warm  mixed 
ocean  surface  layer  of  summer  is  cooled  until  it  reaches  its 
freezing  point  and  ice  crystal  formation  commences.  Ice 
initially  forms  around  the  boundaries  of  the  polar  ice  pack 
and  over  the  shallow  protected  waters  of  high-latitude 
coastlines.  The  ice-covered  areas  continue  to  expand  until 
they  merge,  forming  an  ice-locked  Arctic  Ocean  from  October 
to  June.  The  ice  cover  starts  from  a  summer  minimum  of 
approximately  5.2  x  10*^  km^  then  more  than  doubles  in  areal 
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extent  to  a  maximum  of  11.7  x  10  km  by  the  end  of  the  ice 
season  (CIA,  1978,  p.  12). 

The  ice  cover  grows  continuously  until  spring.  Once  the 
days  start  to  lengthen,  the  amount  of  solar  radiation 
increases  until  the  surface  energy  budget  is  again  positive. 
As  noted  earlier,  snow  cover  will  reflect  approximately  80 
percent  of  the  sun's  incident  radiation,  which  slows  the 
initial  heating  stages.  However,  once  the  air  temperature 
reaches  the  snow's  melting  point,  the  albedo  rapidly  drops  to 
50  percent  and  lower.  This  results  in  a  period  of  rapid  snow 
melt.  Continued  heating  initiates  melting  of  the  underlying 
ice,  causing  cracks  and  flaws  to  develop.  The  surface  melt 
water  drains  through  the  cracks,  further  eroding  them  until 
the  ice  breaks  up  into  floes.  Eventually,  the  ocean  surface 
layer  heats  up,  the  floes  gradually  dissolve,  and  the  cycle 
is  complete  (Langleben,  1972). 

Extremes  of  ice  cover  in  summer  and  winter  for  the  period 
1953-1977  were  shown  in  Figure  2.1.  Obviously,  the  sea  ice 
cover  experiences  a  large  degree  of  variability  both 
seasonally  and  interannually .  Ice  cover  is  something  of  a 
misnomer  as  the  sea  surface  is  rarely  covered  by  an  unbroken 
expanse  of  ice.  Even  in  the  very  thick  ice  regions  of  the 
winter  polar  pack,  infrared  measurements  indicate  that  up  to 
ten  percent  of  the  area  of  the  Arctic  ocean  consists  of  open 
water  or  thin  ice  from  recently  refrozen  leads  (Sater  et  al., 
1971,  p.  41) . 
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Multi-year  ice  (ice  which  has  survived  a  summer's  melting) 
comprises  the  majority  of  the  polar  pack,  which  averages  3  m 
in  depth  (Bourke  and  Garret,  1987) .  First-year  ice  rarely 
grows  to  a  thickness  greater  than  2  m.  However,  the  depth  of 
the  ice  in  any  location  is  largely  dependent  on  the  external 
forces  of  wind  and  currents.  These  cause  the  ice  to  converge 
and  diverge.  When  a  region  of  ice  converges,  it  buckles, 
folds,  and  overlaps,  forming  a  rugged  terrain  and  areas  of 
considerably  thick  ice.  For  example,  the  Beaufort  Gyre's 
anticyclonic  flow  causes  ice  to  converge  along  the  north  coast 
of  Ellesmere  Island  and  Greenland.  The  number  of  ridges  in 
this  region  is  well  above  the  average  for  the  Arctic  pack 
(Weeks  et  al.,  1971).  The  mean  ice  thickness  here  is  of  the 
order  of  6  m  to  8  m  (Bourke  and  Garrett,  1987)  .  In  contrast, 
the  ice  pack  east  of  Spitsbergen  (Svalbard)  is  not  confined 
by  land  and  is  free  to  diverge.  Here,  the  average  ice 
thickness  is  significantly  less,  averaging  approximately  2  m. 

The  thickness,  distribution  and  concentration  of  sea-ice 
are  influenced  by  both  dynamic  and  thermodynamic  factors.  The 
ice  portion  of  an  accurate  ice-ocean  linked  model  needs  to 
include  as  much  of  the  physics  of  these  mechanisms  as 
possible. 
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III.   COUPLED  NUMERICAL  MODEL 

The  coupled  model  used  in  this  study  contains  an  ice  model 
with  both  dynamic  and  thermodynamic  forcing  linked  to  a  fully 
prognostic  primitive  equation  ocean  model.  The  grid  domain 
includes  the  Norwegian  and  Greenland  Seas,  and  the  entire 
Arctic  basin.  The  geography  is  smoothed  due  to  the 
limitations  of  the  resolution  leaving  only  Spitsbergen  and 
Iceland  as  true  islands.  Bottom  topography  is  included.  The 
model  is  forced  by  wind  stress  and  energy  balances  at  the  sea 
surface  and  water  mass  balance  at  open  boundaries.  The  ocean 
model  is  coupled  with  the  ice  model  by  momentum,  heat,  and 
salt  exchanges  through  a  common  ocean  mixed  layer.  Much  of 
this  chapter  follows  the  very  clear  descriptions  of  Van 
Ypersele  (1986)  who  applied  a  similar  model  to  the  Antarctic. 

A.   OCEAN  MODEL 

The  physical  state  of  the  ocean  is  characterized  by  seven 
quantities:  pressure,  potential  temperature,  salinity, 
density,  and  the  three  components  of  velocity.  A  three- 
dimensional  ocean  model  must  be  able  to  predict  the  evolution 
of  those  seven  quantities  from  an  initial  state  when  adequate 
boundary  conditions  are  specified. 

Seven  equations  are  required  to  relate  the  seven  variables 
and  close  the  system.   These  equations  are  derived  from  basic 
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principles  including  the  conservation  of  momentum,  energy  and 
mass  as  well  as  Newton's  law  of  gravitation.  An  equation  of 
state  is  used  to  calculate  density  from  temperature,  salinity, 
and  depth.  The  equations  must  be  simplified  yet  still  retain 
the  most  important  terms  in  order  to  solve  them  numerically 
and  produce  reasonable  results.  Following  Semtner  (1974, 
1986b)  the  necessary  simplification  was  provided  by  making  the 
following  approximations  and  assumptions  about  the  system  of 
equations: 

-  The  ocean  is  incompressible.  This  eliminates  acoustic 
waves  from  the  possible  solutions  and  leads  to  a  simpler 
continuity  equation  without  density  terms.  It  is  a 
reasonable  assumption  because  the  volume  changes  due  to 
pressure  are  small.  For  example,  the  change  of  pressure 
corresponding  to  a  change  of  depth  of  1000  m  would  alter 
the  volume  of  a  typical  sample  of  seawater  by  less  than 
0.5%.  (Pond  and  Pickard,  1983). 

-  The  ocean  behaves  as  a  horizontally  isotropic  fluid  with 
constant  kinematic  eddy-viscosity  coefficients  in  the 
vertical  and  horizontal  directions.  This  means  that  an 
error  is  introduced  in  the  mixing  of  momentum  if  the 
isopycnals  are  not  horizontal. 

-  Exchanges  of  heat  and  salt  which  occur  at  sub-grid  scale 
can  be  represented  by  similar  eddy-dif fusivity 
coefficients  in  the  vertical  and  horizontal  directions. 
On  the  molecular  scale,  salt  diffusion  in  water  is  on  the 
order  of  100  times  slower  than  heat  diffusion.  However, 
turbulent  processes  tend  to  dominate  which  permits  use  of 
the  same  eddy-dif fusivity  coefficients  for  both  salt  and 
heat. 

-  The  thin-shell  approximation  is  made.  This  allows  certain 
terms  in  the  horizontal  momentum  equations  to  be  neglected 
because  the  depth  of  the  ocean  is  negligible  compared  to 
the  radius  of  the  earth. 

-  The  hydrostatic  assumption  is  made.  The  differential  of 
pressure  with  respect  to  depth  is  considered  only  a 
function  of  the  product  of  density  and  the  acceleration 
of  gravity.    This  approximation  is  a  result  of  scale 
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analysis  of  the  vertical  momentum  equation  which  reduces 
to  the  form: 


=  -Pg 


(3.1) 


-  The  Boussinesq  approximation  is  made.  This  is  basically 
the  neglect  of  small  density  variation  effects  on  the 
fluid  momentum  and  therefore  on  the  acceleration  for  a 
given  force,  while  maintaining  these  same  effects  on 
weight  in  the  earth's  gravity  field  (buoyancy). 
Consequently,  average  density  over  the  domain  can  be  used 
in  the  horizontal  momentum  equations.  However,  the  actual 
density  values  are  required  to  derive  pressure  from 
Equation  3.1. 

-  The  Coriolis  terms  and  the  viscous  terms  involving  the 
vertical  speed,  w,  in  the  horizontal  momentum  equations 
are  neglected  on  the  basis  of  scale  analysis. 

-  The  rigid  lid  approximation  is  made.  This  is  achieved  by 
specifying  that  the  vertical  velocity,  w,  be  zero  at  the 
upper  surface.  This  process  filters  out  high-speed 
surface  gravity  waves  in  order  to  permit  use  of  a  longer 
timestep. 

The  three  momentuTn  equations  are  derived  from  Newton's 
second  law  of  motion.  The  resultant  equations  in  vector  form 
are  often  called  the  Navier-Stokes  equations.  They  provide 
a  relation  between  the  velocity,  the  pressure  gradient  force, 
the  viscous  stresses,  and  the  body  forces  such  as  gravity  and 
Coriolis  (e.g.,  Tritton,  1977). 

After  the  simplifications  are  made,  as  allowed  by  the 
assumptions  and  approximations  described  above,  the  momentum 
equations  in  spherical  coordinates  become: 
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where  f  is  the  Coriolis  parameter  (f  =  2nsin0[s"^  ] )  ,  fi  is  the 
angular  speed  of  rotation  of  the  earth,  (n  =  7.272  x  10' 
rad/s)  .  A^  and  A^^  are  the  vertical  and  horizontal  eddy 
viscosity  coefficients,  respectively  (0.3  x  10'^  mVs  and  12 
X  10^  mVs)  .  r^  is  the  earth's  radius  (6.37  x  10®  m)  .  L(a) 
is  the  advection  operator  which  in  spherical  coordinates 
becomes : 
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The  horizontal  Laplacian  operator  becomes: 
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(3.6) 


where  a  is  an  arbitrary  quantity. 

Equation  3.4  is  the  hydrostatic  equation  governing  the 
depth  related  variation  of  pressure  in  the  ocean.  It  assumes 
that  a  perfect  balance  of  forces  exists  in  the  vertical  (i.e., 
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vertical  accelerations  are  negligible) .  This  is  valid  for 
large  scale  motions,  but  it  is  not  adequate  to  describe  small- 
scale  convection,  where  horizontal  and  vertical  scales  of 
motion  tend  to  be  similar.  Therefore,  the  model  must  be 
supplemented  by  an  implicit  treatment  of  convection,  which 
will  prevent  the  unrealistic  condition  of  unstable 
stratification.  The  convective  adjustment  scheme  adopted  here 
is  such  that  if  the  condition 

/Xa.S,p)K  >  p(Q,S,p)ui  (3.7) 

is  verified,  then  the  new  values  of  temperature  and  salinity 
are  recomputed  from 

0k  =  0k.i=¥    and   Sk  =  Suits'  (3.8) 

Level  k  +  1  is  taken  to  be  the  common  reference  level,  and  (~) 
indicates  a  weighted  average  over  layers  k  and  k  +  1. 

The  conservation  of  mass,  together  with  the  assumption  of 
incompressibility ,  gives  the  continuity  equation: 


V  ■  V  =  0 
In  spherical  coordinates  this  becomes 


(3.9) 
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The  theory  of  thermal  transfer  in  a  fluid  of  constant  heat 
capacity,  together  with  an  eddy-dif fusivity  closure 
hypothesis,  gives  a  conservation  equation  relating  the 
evolution  of  temperature  to  the  motion  field  and  the  first 
and  second  spatial  derivatives  of  the  potential  temperature: 

^+L(B)=  Kz^+KhVh0  (3.11) 


Kj  and  K^  are  the  vertical  and  horizontal  eddy  diffusivity 
coefficients,  respectively  (0.3  x  10''^  mVs  and  2  x  10^  mVs)  . 
A  similar  equation  is  derived  for  salinity: 


^     <     1  (Q\  T/  d"^^     <     T    v'^Q  (3.12) 


An  equation  of  state  is  used  to  express  the  density  as  a 
function  of  potential  temperature,  salinity,  and  pressure: 

For  the  purposes  of  numerical  modelling,  the  equation  of 
state  is  expressed  as  a  polynomial  approximation  to  the 
results  of  laboratory  experiments.  The  approximation  of 
Eckart  (1958) ,  utilized  by  Semtner  (1974)  will  also  be  used 
for  this  study. 
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B.   ICE  MODEL 

The  sea-ice  portion  of  the  coupled  model  contains  both 
thermodynamics  and  dynamics.  The  thermodynamic  portion  is 
based  on  the  three-layer  thermodynamic  ice  model  proposed  in 
Semtner  (1976a) ,  which  is  a  simplification  of  the 
comprehensive  ice  model  of  Maykut  and  Untersteiner  (1971). 
It  is  forced  from  above  by  the  atmospheric  fluxes  computed 
from  atmospheric  data,  and  from  below  by  oceanic  heat  fluxes 
prescribed  by  the  ocean  model.  The  dynamic  portion  is  a 
simplification  of  the  "viscous-plastic"  rheology  approach  used 
by  Hibler  (1979).  The  simplification  reduces  the 
computational  load  significantly  and  permits  use  of  the  same 
monthly  averaged  atmospheric  data  for  forcing  as  used  by  the 
ocean  model. 

For  the  three-layer  thermodynamic  model,  sea-ice  is 
assumed  to  be  a  horizontally  uniform  slab  of  ice,  represented 
by  two  layers  of  equal  thickness.  A  covering  layer  of  snow 
is  possible.  When  ice  is  present,  the  temperature  in  the 
oceanic  mixed  layer  below  the  ice  cover  is  assumed  to  be  at 
the  freezing  point  of  seawater  (-1.9°C).  The  temperature 
profile  through  the  ice  and  snow  is  governed  by  the  one- 
dimensional  heat  equation: 


l'^-Ui=^-(^  (3.14) 
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where  p  is  the  density  [900  and  330  kg  m'''  for  ice  and  snow, 
respectively],  C  is  the  specific  heat  capacity  [2.09  kJ  kg'^ 
K'^  for  both  ice  and  snow] ,  T  is  the  temperature  [K] ,  and  k 
is  the  thermal  conductivity  [4.068  and  0.31  W  m'^K'^  for  ice 
and  snow,  respectively] . 

At  the  top  and  bottom  surfaces  of  the  cover  (ice  or  snow 
as  applicable) ,  a  balance  of  fluxes  is  forced  to  exist  such 
that  any  heat  excess  or  deficit  that  occurs  at  either  boundary 
is  used  to  melt  or  freeze,  respectively.  Melting  can  occur 
at  either  surface;  however  freezing  can  only  occur  at  the  ice 
bottom.  Therefore  when  a  negative  flux  imbalance  occurs  at 
the  surface,  the  heat  equation  is  applied  to  the  ice  to 
ascertain  the  effect  at  the  ice  bottom  before  the  amount  of 
freezing  is  determined.  The  melting  point  for  snow  is  set  at 
0°C  and  for  ice  -0.1°C.  The  heat  of  fusion  for  snow  is  set  at 
110  MJ/m^  and  for  ice  at  301  MJ/m^.  Diffusive  fluxes  must  be 
equal  at  the  internal  snow-ice  interface. 

By  late  summer,  the  ice  will  have  completely  melted  in 
certain  regions.  A  heat  budget  equation  is  then  applied  to 
the  vertically  isothermal  surface  mixed  layer  and  its 
temperature  is  allowed  to  increase.  Horizontal  mixing  and 
advection  of  the  warmed  oceanic  surface  layer  can  then  spread 
the  heat  under  adjacent  ice  covered  areas  to  speed  the  melting 
process  there.  Once  the  freezing  season  commences,  the  heat 
budget  cools  the  mixed  layer  to  its  freezing  point,  ice  starts 
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to  form  and  the  heat  equation  for  ice/snow  must  again  be 
applied. 

In  order  to  determine  the  final  temperature  of  the  upper 
surface  when  ice  is  present,  several  conditions  are  imposed 
to  control  the  way  the  various  fluxes  are  applied.  If  the 
initial  temperature  at  the  upper  surface,  Tj.^^^^  <  T^^^  ,  the 
upper  surface  flux  balance  equation  is  given  by: 


Qsw(l-Io)  (l-asurr)  +  Qfw  "  ^Iw  "^  ^11  +  ^LE  +  ^coud  =  ° 


(3.15a) 


where  a^^^^^  is  the  surface  albedo  and  Iq  is  the  fraction  of 
solar  radiation  absorbed  in  the  ice  cover.  Iq  is  set  equal 
to  0  for  snow  and  to  0.17  for  bare  ice.  The  remainder  of  the 
terms  are: 

-  Q^y  longwave  heat  flux,  the  arrow  indicates  downward  (+ve) 
or  upward  (-ve) . 

-  Q^   sensible  heat  flux  (+ve  downward) . 

-  Qle   latent  heat  flux. 

-  Qgy   short  wave  solar  flux. 

-  Qcond   flux  conducted  downward  from  the  surface  to  the 
bottom  of  the  ice. 

These  are  described  in  more  detail  in  section  E  (Forcing 

Fields) . 

If  Tg^^^  =  T^^y^    ,    the  upper  surface  flux  balance  equation  is 

given  by: 
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where  hj  and  h^  are  the  thicknesses  of  the  ice  and  snow, 
respectively,  and  Ly,  is  the  latent  heat  of  fusion. 

At  the  bottom  of  the  ice,  the  flux  balance  equation  is: 


coad  ocean  —  ^Mr^lJi' 


Z  =  \\\ 


(3.16) 


-2. 


where  Qocean  ^^  "^^^  oceanic  heat  flux  [W  m'  ]  . 

The  outgoing  longwave  radiation  Q|_y  follows  the  Stefan- 
Boltzmann  law: 


'^^^        hiirf 


(3.17) 


where  e  is  the  surface  emissivity  (non-dimensional,  and 
assumed  equal  to  1)  ,  and  o  is  the  Stef an-Boltzmann  constant 
[=  5.67  X  10"^  W  m'^  K""*]  . 

The  final  temperature  at  the  upper  surface,  T^^^^^,  is 
obtained  by  linearizing  the  blackbody  emission  term  3.17. 

The  conductive  flux  Qcond  inside  the  snow  and  ice  is 
computed  by  assuming  a  linear  vertical  temperature  profile 
between  grid  points,  so  that  at  the  bottom  of  the  ice: 
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where  kj  is  the  ice  thermal  conductivity,  T^yfreez  ^^  ^^^  freezing 
point  of  seawater  [=271.3  K]  ,  Tj  is  the  temperature  of  the 
second  layer  of  ice,  and  h-  is  the  ice  thickness. 

The  storage  of  latent  heat  in  brine  pockets  is  modeled  by 
an  internal  heat  reservoir  as  in  Semtner  (1976a).  The 
reservoir  accumulates  a  fraction  of  the  solar  radiation  that 
penetrates  the  snow-free  ice  [Q  (1-Iq)  (l-a^^g)  ]  .  Total  heat  in 
the  reservoir  is  limited  to  50%  of  the  heat  needed  to  melt  all 
the  ice  at  that  position.  This  stored  heat  adds  inertia  to 
the  simulated  melting  process.  The  effects  of  rapid  changes 
in  thermodynamic  forcing  are  smoothed  and  the  onset  of  the 
melting  season  is  delayed.  A  schematic  of  the  fluxes  and 
processes  incorporated  in  the  thermodynamic  portion  of  the  ice 
model  is  shown  in  Figure  3.1.  The  numerical  scheme  used  to 
compute  the  ice  and  snow  thickness  and  the  temperature  at  the 
different  levels  is  contained  in  Semtner  (1976a)  and  will  not 
be  repeated  here. 

The  ice  dynamics  portion  of  the  model  uses  a  simplified 
Hibler  (1979)  rheology  called  "bulk  viscous"  (Washington  and 
Parkinson,  1986;  Hibler,  1988).  This  approach  was  developed 
for  use  on  seasonal  time  scale  ice  models  with  slowly  varying 
forcing.  Semtner  (1987)  used  this  rheology  and  was  able  to 
produce  annual  ice  cycles  very  similar  to  the  Hibler  and  Bryan 
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Figure  3 . 1 


Schematic  of  fluxes  and  ice  processes  incorporated 
in  the  thermodynamic  portion  of  the  ice  model. 
The  total  heat  flux  provided  by  the  ocean 
circulation  below  the  mixed  layer  is  shown  as  Q„^„,„ 
and  Qcond  ^^  "^^^  amount  of  heat  conducted  through 
the  ice  from  the  surface.  The  heat  reservoir  can 
only  receive  heat  when  the  ice  is  snow  free; 
however,  it  can  provide  heat  whenever  the  net  heat 
balance  at  the  surface  becomes  negative. 
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(1986)  model  which,  as  noted  earlier,  used  the  complete  daily 
variable  forcing  dynamic  approach  of  Hibler  (1979)  .  The 
smoothed  forcing  permits  an  explicit  timestepping  technique 
to  be  used  vice  an  iterative  approach  and  the  stress  tensor 
formulation  is  much  simpler.  These  two  factors  reduce  the 
computational  requirements  by  approximately  66%.  The  dynamics 
portion  of  the  ice  model  predicts  changes  to  the  ice  thickness 
and  concentration  from  advection,  diffusion  (a  numerical 
requirement)  and  compression  effects.  The  sea-ice  velocity, 
V- ,  is  predicted  as  follows: 

DV,-/Dt  =  l/m(-mfVj  +  t^   +  t^  -  mg  Vh  +  F)         (3.18a) 

where  D/Dt  is  the  substantial  time  derivative,  f  is  the 
Coriolis  parameter  and  F  is  the  internal  ice  stress.  t^  and 
r^  are  the  wind  and  ocean  current  stresses,  respectively  and 
Vh  is  the  sea  surface  slope.  m  is  the  mass  and  g  is  the 
acceleration  of  gravity. 

The  bulk  viscous  rheology  causes  the  ice  to  resist 
compression  in  a  plastic  fashion  but  allows  unimpeded 
divergence.  In  the  Hibler  (1979)  notation  the  stress  tensor 
is  reduced  to: 


O  : 


..   =   r    e-    6..  (3  .  18b) 


where  the  strain  rate  tensor  e^j  is 


^i  =  5[li-|^j  <-"=' 
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The  bulk  viscosity  C  is  set  to  zero  for  divergence  (e^j>0) . 
For  convergence,  f  =  P*/2e,-j,  where  P*  =  h,-  x  10*  [N/m^] .  The 
bulk  viscosity  is  not  allowed  to  exceed  2.5  x  10®  x  P*  [sec]. 

C.   GRID 

The  ocean  model  uses  a  spherical  coordinate  system  with 
the  equator  of  the  system  passing  through  the  geographic  North 
Pole  and  the  prime  meridian  situated  at  40  degrees  East 
longitude.  Shifting  the  poles  of  the  spherical  coordinate 
system  90  degrees  to  the  geographic  equator,  removes  the 
singularities  of  the  system  and  provides  more  uniform 
horizontal  resolution  within  the  Arctic  region.  The  resulting 
horizontal  grid  interval  is  approximately  110  kilometers.  The 
position  of  a  grid  point  is  defined  by  latitude,  longitude  and 
ocean  depth. 

The  ice  portion  of  the  linked  model  is  run  on  a  cartesian 
grid  with  110  km  spacing.  The  two  grids  coincide  exactly  at 
the  pole  but  are  slightly  offset  at  the  southern  boundaries. 
This  introduces  a  small  error  in  the  ice  mass  balance  and  salt 
fluxes  into  the  ocean  in  these  areas.  Hibler  and  Bryan  (1987) 
determined  that  the  maximum  errors  introduced  into  their  model 
using  a  similar  scheme  were  only  10%.  A  10%  error  at  the 
southern  fringes  of  the  grid  was  not  considered  critical  to 
the  results  of  this  work. 

The  gridded  domain  used  in  this  study  was  shown  in  Figure 
1.1.   The  grid  points  are  staggered  on  a  lattice  of  type  B 
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(Mesinger  and  Arakawa,  1976) ,  where  the  temperature,  salinity, 
and  stream  function  points  are  located  at  the  center  of  a  grid 
box.  The  horizontal  components  u  and  v  of  the  velocity  are 
at  the  middle  of  the  vertical  edges.  The  vertical  component 
w  is  computed  at  the  center  of  the  horizontal  faces  for  the 
determination  of  heat  and  salt  vertical  advection,  and  at  the 
box  corners  for  the  calculation  of  u  and  v  (Figure  3.2). 

A  total  of  13  levels  are  used  in  the  vertical  in  order  to 
resolve  the  bathymetric  variations  shown  in  Figure  3.3.  The 
uppermost  level  is  treated  as  an  isentropic  mixed  layer  of  30 
m.  The  levels  are  thinnest  near  the  surface  to  allow  five  of 
these  levels  above  the  continental  shelf,  which  is  at 
approximately  250  m  depth.  Three  more  levels  are  defined  down 
to  approximately  1100  m  depth  for  representation  of  the 
Atlantic  layer.  Table  1  lists  the  13  layer  thicknesses  and 
the  depth  at  the  bottom  of  each  layer. 

D.   BOUNDARY  AND  INITIAL  CONDITIONS 

The  rigid  lid  approximation  requires  that  w=0  at  the  upper 
boundary  (z=0).  The  wind  stress  r [Pa]  is  specified  at  the 
surface  by 

r{^..(t>)  =  PoAzj^iu^y)    at  z  =  0  (3.19) 

The  small  amount  of  water  mass  which  is  transferred  at 
the  ocean  surface  is  neglected.  This  includes  precipitation 
minus  evaporation,  ice  melt  and  minor  river  runoff.   However, 
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Figure  3.2  The  distribution  of  levels  and  grid  points. 
Labels  to  the  right  and  left  of  the  grid  points 
refer  to  ocean  and  sea-ice  variables, 
respectively . 
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Figure  3 . 3 


Arctic  bathymetry  plotted  with  a 
interval.   (froin  Semtner,  1987). 


500  in  contour 
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TABLE  1 
VERTICAL  LAYER  THICKNESSES  AND  CORRESPONDING  DEPTHS 


LAYER 

THICKNESS  (M) 

DEPTH  (M) 

1 

30 

30 

2 

20 

50 

3 

40 

90 

4 

60 

150 

5 

100 

250 

6 

183 

433 

7 

299 

732 

8 

391 

1123 

9 

487 

1610 

10 

529 

2139 

11 

563 

2702 

12 

585 

3287 

13 

589 

3876 

the  resultant  changes  in  salinity  from  these  processes  are 
accounted  for  by  applying  prescribed  negative  salt  fluxes  at 
the  surface.  The  prescribed  salt  fluxes,  as  well  as  the  salt 
flux  introduced  during  simulated  freezing,  are  accounted  for 
by  the  salt  flux  density  Q.[K  m'^  s'^]: 


dS 


Q^  =  PwKz-^  =^  -pwSSVV 


(3.20) 


where  SSW  [m/s]  is  the  sum  of  all  fluxes  of  fresh  water  at 
the  surface. 

A  heat  flux  density  Q^[Wm'^]  specified  at  the  surface  (z=0) 
affects  the  temperature  field  in  the  same  way: 


Qx  =  Kz§  -  SSH 


(3.21) 
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where  SSH  [W/m^]   is  the  net  energy  flux  density  at  the 
surface. 

Vertical  motion  is  forced  at  the  ocean  bottom  as  the 
currents  respond  to  the  vertical  contours  of  the  bottom  slope. 
Bottom  friction  is  determined  using  a  quadratic  drag  law  in 
the  same  manner  as  Weatherly  (1972) .  Vertical  advective 
fluxes  of  heat,  salt  and  momentum  at  the  bottom  are  taken  to 
be  zero.  The  vertical  diffusive  fluxes  of  heat  and  salt  at 
the  bottom  are  also  zero. 

PoKz-^  (T,S)  =  0      atz  =  -H(A,0)  (3.22) 

The  horizontal  boundaries  coincide  (to  the  limits  of  the 
resolution)  with  the  geographic  coastlines  and  islands  of  the 
Arctic  Ocean  and  adjacent  seas.  A  no-slip  condition,  i.e., 
(u,v)=0,  is  imposed  at  lateral  walls,  with  no  heat  or  salt 
fluxes  across  them. 

Each  model  run  simulates  the  ocean  and  ice  cover  for  one 
year,  producing  output  fields  for  the  end  of  each  month.  The 
decade  from  1971-1980  was  simulated  by  conducting  ten 
successive  runs  using  the  final  (end  December)  ocean  and  ice 
condition  from  the  previous  year's  simulation,  as  the  start 
condition  for  the  current  simulation.  Since  the  analysis 
begins  with  1971,  December  1970  would  be  used  as  the  very 
first  start  condition.  However,  because  the  model  requires 
some  time  to  achieve  a  steady  state  after  each  modification, 
a  three  year  run-up  was  used.   That  is,  the  runs  were  started 
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at  1968,  giving  the  model  a  chance  to  stabilize  before  the 
data  from  1971-1980  were  produced  and  subsequently  extracted 
for  analysis. 

E.   FORCING  FIELDS 

This  section  describes  the  various  forcing  fields  and  data 
used  to  drive  the  coupled  model.  Several  formulas  are 
presented  which  contain  empirical  constants.  While  the  values 
used  for  these  constants  represent  current  estimates,  their 
true  values  remain  uncertain. 

Monthly  averages  of  observed  atmospheric  surface  pressure, 
specific  humidity,  long  and  short  wave  radiation  and  air 
temperature  were  provided  by  Dr.  J.  Walsh.  The  pressure 
fields  were  obtained  from  the  NCAR  daily  analysis  of  surface 
pressure.  This  analysis  is  based  on  daily  reports  from 
stations  on  the  periphery  of  the  Arctic  Basin  and  several  ice 
stations.  Recent  comparisons  of  these  analyses  with  Arctic 
Buoy  data  show  that  the  average  error  is  limited  to  1-2  mb 
(Walsh,  private  communication) .  Temperature  fields  were  also 
obtained  from  periphery  meteorological  stations  and  ice 
stations  however  they  were  compiled  as  monthly  averages.  The 
temperatures  were  assumed  to  be  10  m  values.  Specific 
humidity  at  10  m,  q^Q,  was  calculated  using  the  prescribed 
temperature  fields  and  an  assumed  relative  humidity  of  90%. 
Specific  humidity  at  the  surface,  q^,  was  calculated  using  the 
surface  temperature  produced  by  the  model   and  assuming 
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saturation.  The  long  and  shortwave  radiation  fields  were 
determined  in  a  similar  manner  to  Parkinson  and  Washington 
(1979)  .  Radiation  was  calculated  for  each  hour  as  a  function 
of  prescribed  zenith  angle  and  cloud  cover.  The  cloud  cover 
varied  monthly  over  an  annual  cycle  but  was  constant  across 
the  domain. 

Since  the  prescribed  cloud  cover  was  probably  the  weakest 
of  the  assumptions  made,  Walsh  has  noted  (private 
communication)  that  the  long  and  shortwave  radiation  fields 
probably  contain  the  most  error  of  all  the  atmospheric  forcing 
data.  The  specific  humidities  are  limited  by  the  10  m  height 
and  90%  relative  humidity  assumptions,  but  are  still 
considered  reasonable.  The  pressure  and  temperature  fields 
appear  quite  accurate  as  confirmed  by  monitoring  stations. 

These  were  the  same  forcing  data  used  in  Walsh  et  al., 
(1985)  .  That  work,  as  well  as  that  of  Hibler  and  Walsh 
(1982),  showed  that  inclusion  of  the  dynamic  influences  from 
interannually  varying  atmospheric  forcing  provided  improved 
simulations  of  the  interannual  variations  in  the  sea-ice 
cover.  Time  interpolation  on  a  daily  basis  was  conducted 
using  a  cubic  spline  method  similar  to  Maykut  and  Untersteiner 
(1971).  The  smoothed  daily  forcing,  as  opposed  to  the 
observed  daily  forcing,  was  used  in  order  to  employ  a 
simplified  Hibler  ice  dynamics  parameterization  as  used  in 
Semtner  (1987).  This  reduced  the  computational  load  of  the 
ice  portion  of  the  model  considerably. 
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wind  stress  on  the  ice  and  open  ocean  was  computed  by 
conventional  aerodynamic  bulk  formula  methods  in  the  same 
manner  as  Semtner  (1987) : 

r^  =  Cd/?a||Vg||Vg[B]  (3.23) 

where  C^  is  the  drag  coefficient  (1.75  x  10'^).  This  value  is 
larger  than  the  value  used  by  Hibler  and  Bryan  (1987)  to 
compensate  for  the  use  of  time-averaged  wind.  B  is  a  rotation 
matrix  which  shifts  the  forcing  vector  25  degrees  to  the  left 
of  the  surface  geostrophic  wind  when  the  stress  acts  on  the 
ice.  When  the  wind  stress  is  applied  to  the  open  ocean,  B  was 
set  to  one  (no  deflection)  .  The  wind  stress  on  the  ice- 
covered  ocean  was  not  computed  explicitly  but  was  computed 
interactively  through  the  ice  (See  Section  G.  Ice-Ocean  Model 
Coupling. ) 

Mass  inflow  through  the  Faroe-Shetland  Channel  and  the 
Bering  Strait  and  matched  outflow  through  the  Denmark  Strait 
and  the  Canadian  Archipelago  were  prescribed.  These  values 
are  the  same  as  in  Semtner  (1987).  The  mass  and  T,S  values 
of  the  inflow,  and  the  mass  of  the  outflow  are  invariant  in 
time,  but  the  T,S  properties  of  the  outflow  are  dependent  on 
the  simulated  fields  produced  by  the  model.  The  inflow  of 
fresh  water  from  major  rivers  around  the  margins  of  the  Arctic 
Basin  is  specified  on  the  basis  of  monthly  values  tabulated 
by  Cattle  (1985).  Cubic  splines  are  used  to  interpolate  to 
daily  values. 
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The  flux  of  sensible  heat,  Q„  (positive  downward) ,  is 
computed  by  the  standard  aerodynamic  bulk  formula  method 
(e.g. ,  Gill,  1982) . 

QH  =  PaCpCj^||Vg||(Ta-Tw)  (3.24) 

where  Cp  is  the  specific  heat  of  air  (1004  J/kg/K)  for  dry 
air  (Huschke,  1959),  C„  is  a  dimensionless  transfer 
coefficient  for  sensible  heat  called  the  Stanton  Number  (1.75 
X  10'  (Maykut,  1978) ,  and  T^  is  the  water  surface  temperature. 
The  flux  of  latent  heat,  Q^^,  is  computed  in  a  similar 
fashion: 


Q^E  "=  /'-^^^E^E"^'sll(^,Q-n.surf)  (3.25) 

where  I^  is  the  latent  heat  of  vaporization  or  sublimation 
(2.58  X  10*^  J/kg)  .  C^  is  a  dimensionless  transfer  coefficient 
called  the  Dalton  Number  (1.75  x  10'^  (Maykut,  1978)),  and  q^^ 
and  qsurf  ^^^  "^^^  specific  humidities  at  10  m  and  at  the 
surface,  respectively. 


F.   SOLUTION 

A  comprehensive  description  of  the  solution  methods  for 
the  ocean  model  is  contained  in  Semtner  (1986b).  The 
following  brief  summary  provides  a  general  overview  of  the 
procedures . 

The  numerical  solution  of  Equations  3.2,  3.3,  3.4,  3.7, 
3.8,  3.9,  3.11,  3.12  and  3.13  is  conducted  in  several  steps. 
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The  horizontal  velocity  is  decomposed  into  a  baroclinic 
(effect  of  vertical  shear)  and  a  barotropic  (vertically 
averaged)  component: 


(u,v)  =  (u',v')  +  (u,v) 


(3.26) 


The  vertically  averaged  velocity  components  (u,v)  are  written 
in  terms  of  a  stream  function  0.  The  vertically  integrated 
flow  is  non-divergent  which  guarantees  the  existence  of  0. 
The  stream  function  represents  the  integrated  mass  transport: 


,      1  di) 


-    1  r' 

-_  1  r°  1     1   di!) 


(3.27) 


A  prediction  equation  for  the  stream  function,  independent  of 
the  surface  pressure  is  then  derived;  first  by  vertically 
averaging  and  then  by  taking  the  curl  of  the  horizontal 
momentum  Equations  3.2  and  3.3: 
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where  G*  and  G^  represent  all  the  nonlinear  and  viscous  terms 
in  Equations  3.2  and  3.3. 

Equation  3.28  essentially  relates  the  time  derivative  of 
the  Laplacian  of  the  stream  function  to  the  curl  of  the 
vertically  averaged  sum  of  all  the  forcing  terms  in  the 
Navier-Stokes  equations.  The  solution  requires  inversion  of 
a  second-order  differential  operator  in  a  closed  or  multiply 
connected  domain.  The  stream  function  0  must  be  a  constant 
along  coastlines  to  avoid  mass  transport  across  the 
boundaries.  0  is  held  constant  in  time  along  the  continent 
boundaries.  On  the  islands,  ^  varies  in  response  to  the 
changing  circulation  in  the  same  manner  as  Takano  (1974),  but 
with  some  modification  to  allow  for  variable  bottom 
topography.  Essentially  this  involves  integrating  the  curl 
of  the  vertically  averaged  momentum  equations  around  each 
island.  The  condition  is  then  applied  so  that  the  line 
integral  of  the  surface  pressure  gradient,  VP^  ,  around  each 
island  is  zero  which  produces  an  equation  to  predict  the 
change  of  0  on  each  island. 

A  prediction  equation  for  the  vertical  shear  of  velocity 
is  found  by  differentiating  the  equation  of  horizontal  motion 
with  respect  to  depth  z,  and  applying  the  hydrostatic  relation 
to  eliminate  the  atmospheric  pressure  P^.  The  baroclinic  and 
barotropic  components  of  the  velocity  are  then  added  to  obtain 
the  total  u  and  v. 
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The  prediction  equations  for  temperature  and  salinity 
(3.11)  and  (3.12)  are  much  easier  to  solve.  A  "leapfrog" 
timestep  is  used  for  advection  and  a  "forward"  timestep  is 
used  for  diffusion  to  maintain  numerical  stability. 
Additionally,  a  forward  timestep  is  applied  every  10  steps 
for  the  advection  calculations.  This  was  an  efficient  method 
to  suppress  the  computational  mode  of  the  leapfrog 
timestepping  technique. 

Two  different  timesteps  were  used  based  on  the  different 
rates  of  change  of  temperature,  salinity,  and  velocity  in  the 
ocean  model.  Temperature  and  salinity  calculations  employed 
6-hour  timesteps  while  the  faster  varying  velocity  field  used 
a  shorter  timestep  of  22.5  minutes.  This  reduced  the  total 
number  of  calculations  required,  thereby  speeding  up  the 
numerical  integration.  The  same  technique  was  used 
successfully  by  Bryan  (1984).  Different  timesteps  were  also 
used  in  the  ice  model.  Thermodynamic  timesteps  were  15 
minutes  and  dynamic  timesteps  were  15  seconds  to  maintain 
numerical  stability  as  in  Semtner  (1987)  . 

G.   ICE-OCEAN  MODEL  COUPLING 

The  ocean  and  sea-ice  portions  of  the  linked  model  are 
coupled  through  momentum,  heat  and  salt  exchanges.  The 
momentum  transfer  involved  with  the  coupling  occurs  at  the 
ice-ocean  interface.  The  ice  to  water  momentum  transfer  is 
modeled  by  a  linear  drag  formula: 
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Av  -  Az^p^  =  Cio(Vi  -  Vw)  (3.29) 

where  Vj  is  the  ice  velocity,  V^^  is  the  ocean  velocity  in  the 

first  (mixed)  layer  and  the  linear  drag  coefficient  Cj^  is 

assigned  a  value  of  27.5  x  lO'"^  following  Semtner  (1987).   If 

both  ice  and  open  water  are  present  in  a  grid  space,  the  wind 

stress  is  applied  to  the  open  water  with  Equation  3.19  and  the 

ice  to  ocean  stress  is  computed  by  Equation  3.29.  The  average 

stress  is  then  computed  by  weighting  the  two  contributions 

appropriately  by  the  ice  concentration. 

The  oceanic  heat  flux  into  the  mixed  layer  is  provided  by 

the   ocean  model.     The   computation  of  the  mixed   layer 

temperature  at  each  time  step  and  grid  point  is  made  in  three 

parts : 

-  The  first  computation  isolates  the  first  layer 
thermodynamically  and  dynamically  from  the  rest  of  the 
ocean.  The  temperature  of  the  layer  T^^ ,  is  determined 
simply  from  the  mixed  layer/atmospheric  heat  balance.  If 
T^^  reaches  the  freezing  point  of  seawater,  sea-ice  is 
formed.  The  change  in  temperature  from  the  preceding  time 
step  is  equivalent  to  a  heat  flux  density  Q^^  (W/m^)  at  the 
surface : 


Q^l  =  PwCwhmi%^  (3.30) 

where  h^^  is  the  depth  of  the  mixed  layer  (30  m) . 

The  second  step  uses  the  Q^^  calculated  in  Equation  3.30 
as  an  upper  surface  thermodynamic  forcing  for  the 
computation  of  vertical  heat  diffusion  in  the  ocean.  The 
ocean  distributes  this  input  of  energy  by  horizontal 
diffusion,  advection  and  convection.  A  new  temperature 
Th,^.^;.  is  then  obtained  for  the  mixed  layer. 

dynamic  -' 
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-  The  final  step  ensures  that  energy  is  conserved.  If  sea- 
ice  is  present  at  the  grid  point  being  worked  with,  and 
if  Tjjy^g^j^  is  higher  than  the  freezing  point  of  seawater, 
T^^n^.^, /  then  enough  ice  is  melted  to  decrease  T^^„„-  back  to 
the  freezing  temperature.  In  the  rare  cases  where  T^ynamic 
is  actually  lower  than  the  freezing  temperature,  ice  will 
accrete  to  the  bottom  of  the  existing  ice  until  the  latent 
heat  of  freezing  provides  enough  heat  to  return  the  layer 
to  the  freezing  temperature. 

Salt  exchanges  are  computed  by  prescribing  the  salt  flux 
based  on  the  amount  of  sea-ice  melting  and  accretion, 
precipitation,  evaporation  and  glacial  melt.  When  no  sea-ice 
is  present,  the  mixed  layer  salinity  is  computed  from  Equation 
3.12  and  3.20  by  specifying  the  fresh-water  flux  at  the 
surface,  SSW. 

When  sea-ice  is  present,  ice  accretion  and  melting  produce 
positive  and  negative  salt  fluxes,  respectively.  It  is 
assumed  that  upon  freezing,  70%  of  the  initial  salt  content 
of  the  seawater  is  rejected  (Semtner,  1987).  Therefore  the 
salt  flux  density  Q^^^^  (kg  s'^m'^)  at  the  upper  surface  (z=0)  is 
given  by: 


Qsalt  =  Kz§  [ps  ^JT^^^ll'c^^  +  0.70/^i^^^]  -/5w[(P-E)(l-Ai)  +  R)    (3.31) 

where  A^  is  the  area  of  the  ice  cover;  a^^^  =  1  if  snow  is 
melting  and  0  otherwise.  P  is  the  precipitation  rate,  E  is 
the  evaporation  rate  and  R  is  the  amount  of  glacial  melt. 

Table  2  is  a  summary  of  the  constants  used  in  the  model. 
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TABLE    2 
CONSTANTS 


Time   steps 


Arakawa  type  B  grid 


ocean  T  and  S 
ocean  velocity 
ice  T 
ice  velocity 


acceleration  of  gravity 

earth's  speed  of  rotation 

earth's  radius 

vertical  eddy  viscosity  (u,v) 

horizontal  eddy  viscosity  (u,v) 

vertical  eddy  diffusivity  (T,S) 

horizontal  eddy  diffusivity  (T,S,ice) 
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IV.   INTERACTIVE  OCEAN  EXPERIMENT 

This  experiment  was  conducted  to  determine  the  importance 
of  incorporating  a  fully  interactive  prognostic  ocean  model 
into  a  linked  ice-ocean  numerical  model.  Many  numerical 
models  of  sea  ice  have  been  developed  which  do  not  include  a 
prognostic  ocean  component.  This  experiment  was  designed  to 
test  the  validity  of  that  omission.  The  importance  of 
including  a  prognostic  ocean  was  judged  with  respect  to  the 
accuracy  of  the  simulation  of  the  interannual  and  annual 
variations  of  sea-ice  cover  in  the  Arctic.  Evaluation  was 
based  on  comparisons  between  observed  fields  of  ice 
concentration  and  simulated  ice  concentration  fields  produced 
by  the  model. 

A.   EXPERIMENTAL  SETUP 

The  linked  ice-ocean  model  was  initially  run  for  ten  years 
with  a  fully  interactive  ocean  and  monthly  varying  atmospheric 
forcing.  All  oceanic  variables  were  saved  on  a  regular  basis 
from  this  first  run,  in  order  to  examine  their  variability  and 
to  calculate  a  ten-year  average  annual  cycle.  The  second  run 
used  the  same  atmospheric  forcing,  integration  period  and  ice 
portion  as  in  the  first  run.  However  the  interactive, 
prognostic  ocean  portion  was  replaced  by  prescribing  the 
average  annual  cycle  from  the  first  run.   The  outputs  from 
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both  runs  were  analyzed  to  determine  each  model's  skill  in 

reproducing   the   observed   annual   cycle   and   interannual 

variations  of  ice  concentration.  Modifications  to  the  Semtner 

(1987)  model  included: 

-  Run  1 

Changing  the  atmospheric  forcing  from  a  30-year 
mean  annual  cycle  to  an  interannually  varying 
cycle  of  observed  monthly  averages. 

Defining  the  initial  conditions  for  each  modelled 
year  as  the  end  condition  for  the  previous  year 
vice  averaged  or  constant  fields.  This  permitted 
consecutive  integration  of  as  many  years  as 
desired. 

-  Storing  and  then  averaging  simulated  oceanic  heat 
flux  and  layer  2  ocean  velocities  at  every  time 
step  for  ten  years.  The  oceanic  heat  flux  saved 
for  each  grid  box  was  the  summation  of  all  the 
heat  the  ocean  model  provided  to  the  mixed  layer 
in  that  box.  Layer  2  velocities  (corresponding 
to  a  depth  of  4  0  m)  were  used  because  they 
simulated  the  geostrophic  currents  quite  well. 
The  current  velocities  in  the  mixed  layer  were 
unrealistic  because  all  the  induced  Ekman  flow  was 
confined  within  that  layer.  A  new  data  file 
containing  the  10-year  average  annual  cycle  of 
these  variables  was  then  created. 


-  Run  2 


Incorporating  the  10-year  average  annual  cycle  of 
oceanic  heat  flux  and  layer  2  velocities  into  the 
model's  spline  smoothing  routine.  This  permitted 
the  continued  use  of  the  modified  Hibler  (1979) 
ice  dynamics. 

Replacing  the  ocean  model  with  a  parameterization 
of  the  average  oceanic  forcing  (thermodynamic  and 
dynamic)  .  This  forcing  was  applied  to  the  base 
of  the  mixed  layer. 
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B.   RESULTS 

The  average  computing  time  for  a  single  simulation  year 
using  the  fully  prognostic,  interactive  ocean  in  run  1  was  575 
seconds  on  the  CRAY  XMP.  The  average  time  for  a  single  year 
in  run  2  with  the  prescribed  ocean  was  275  seconds.  For 
labeling  and  descriptive  purposes,  the  observed  data  were 
designated  as  A  data.  Run  1  output  as  B  data  and  Run  2  output 
as  C  data. 

The  simulated  ice  concentration  contours  for  all  120 
months  were  produced  for  B  and  C  in  the  same  manner  as  for  the 
observed  fields.  As  representative  examples.  Figures  4.1  - 
4.6  shov;  the  ice  concentration  contours  for  months  76  and  93 
(Apr  77  and  Sep  78)  for  all  three  data  sets  (A,  B  and  C)  . 
These  contours  are  typical  for  minimum  and  maximum  ice  extent 
periods . 

Layer  2  velocities  (Figures  4.7  -  4.10)  and  contours  of 
oceanic  heat  flux  into  the  mixed  layer  (Figures  4.11  -  4.14) 
for  each  month  for  B  and  the  corresponding  average  fields  for 
C  are  also  shown.  Heat  units  are  degrees  C  per  second, 
applied  to  the  volume  of  the  mixed  layer  within  each  grid  box 
(3.63  X  10^^  m^)  . 

Ice  thickness  was  another  useful  predicted  variable.  The 
ice  thickness  distribution  provided  a  good  indication  of  how 
the  ice  dynamics  in  combination  with  the  ice  rheology  modified 
the  ice  cover  within  the  high  ice  concentration  regions. 
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Figure  4.1  Observed  ice  concentration  (A)  in  tenths  for  month 
76  (April,  1977).  Contour  line  labelled  -1.5  is 
coastline . 


74 


I  I  I  I  I  I  I  I 


1  I  i  I  i  TTl  i  NJ  I  i  I  I  i  1  I  i  I  i  i  I  I  1  Tl  I  I  i  I  I  I 


7 


^.. 


-7.t« 


\ 


-i.e* 


/' 


1  1  I  !  I  I  I  I  !  !  I  '  !  I  1  1  I  !  I  I  I  I  !  I  !  !  I  !  I  !  !  I  r  I  I  I  I  1  !  '  1 


Figure  4.2  Simulated  ice  concentration  using  prognostic  ocean 
mcpdel  (B)  for  month  76  (April,  1977).  Contour 
line  labelled  -1.6  is  coastline.   Contours  are  in 

tenths. 


75 


"T  I  I  I  i  I  I  i  I  1  I  I  I  I  I  I  r?  i  fs~i  I  I  I  I  I  r~i  1  !  \   I  I  I  I  I  I  I  1  I  i  I  I 


Figure  4.3  Simulated  ice  concentration  using  10-year  mean 
ocean  data  (C)  for  month  76  (April,  1977). 
Contour  line  labelled  -1.6  is  coastline.  Contours 

are  in  tenths. 
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Figure  4.4   Observed  ice  concentration  (A)  in  tenths  for  month 
93  (September,  1979)  . 
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Figure  4.5  Simulated  ice  concentration  using  prognostic  ocean 
model  (B)  for  month  93  (September,  1979) . 
Contours  are  in  tenths. 
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Figure  4.6  Simulated  ice  concentration  using  10-year  mean 
ocean  (C)  for  month  93  (September,  1979)  . 
Contours  are  in  tenths. 
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Figure  4  . 7 


10-year  mean,   level   two   (40  m  depth)   ocean 
currents  for  April.  Maxiinum  vector  0.246  m/s. 
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Figure  4.8  Simulated  level  tv.'o 
for  month  76  (April, 
m/s . 


(40  in  depth)  ocean  currents 
1977).   Maximum  vector  0.162 
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Figure    4 . 9 


10-year  nean  level  tv/o    (4  0  m  depth)    ocean  currents 
for   September.       Kaxiraum   vector    0.152    m/s. 

82 


\ 


/  ^r — , 

^   ^  *    -    \   /   y    <     \ 


/  I    '    '   >>  >^  "  '   '  y   / 


1 


'/: 


Figure  4.10  Simulated  level  two  (40  m  depth)  ocean  currents 
for  month  93  (September,  1979) .  Maximium  vector 
0.344  m/ s . 
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Figure  4.11  10-year  mean  ocean  heat  flux  through  the  30  m 
mixed  layer  for  April.  Units  are  degrees 
C/sec/cm^  scaled  by  1  x  10®.  Conversion  factor  to 
W/m^  is  1.25  X  10^.  L  and  H  indicate  relative  lows 

and  highs  respectively. 
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Figure  4.12  Simulated  ocean  heat  flux  through  the  30  m  mixed 


layer  for  month  76.   Units  are  degrees  C/sec/cm 

scaled  by  1  ■•  " '^^    -" --  ------  --  - --^ 

1.25  X  10^ 


X  10".   Conversion  factor  to  W/m  is 
L  and  H  indicate  relative  lows  and 
highs  .respectively. 
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Figure  4.13  10-year  mean  ocean  heat  flux  through  the  30  m 
mixed  layer  for  Sep.   Units  are  degrees  C/sec/cm 

scaled  by  1  x  10^.  Conversion  factor  to  W/m  is 
1.25  x  10^.  L  and  H  indicate  relative  lows  and 
highs  respectively. 
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Figure  4.14  Simulated  ocean  heat  flux  through  the  30  m  mixed 
layer  for  month  93.  Units  are  degrees  C/sec/cm^ 
scaled  by  1  x  10^.  Conversion  factor  to  W/m  is 
1.25  X  10^.  L  and  H  indicate  relative  lows  and 
highs  respectively. 
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Representative  ice  thickness  contours  for  B  and  C  in  months 
76  and  93  are  shown  in  Figures  4.15  -  4.18. 

In  order  to  gain  a  more  quantitative  appreciation  of  the 
differences  between  the  various  outputs  (A,  B  and  C)  ,  a  series 
of  comparative  graphs  were  produced.  The  total  monthly  ice 
area  for  each  of  the  four  regions  was  determined  by 
multiplying  the  ice  concentration  at  each  grid  point  by  the 
grid  box  area  and  summing  over  the  region  (Figures  4.19  - 
4.22).  The  observed  ice  concentration  (A)  was  subtracted  from 
the  B  and  C  fields  to  show  the  amount  each  simulated  field 
differed  from  the  observed  (Figures  4.23  -  4.26).  The 
absolute  sum  of  the  areal  differences  for  B  and  C  in  each 
region  was  then  calculated.  A  measure  of  the  improvement  in 
simulating  the  annual  cycle  of  total  ice  area  in  each  region 
was  then  possible.  The  absolute  sum  of  B-A  was  consistently 
less  than  the  absolute  sum  of  C-A  indicating  that  the  B  data 
were  closer  to  the  observed  data  and  on  average  more  accurate 
than  the  C  data.  The  results  from  these  calculations  are 
shown  in  Table  3. 

Figures  4.19  -  4,22  indicated  that  both  B  and  C  had 
obvious  biases  in  total  ice  area.  In  order  to  examine  the 
interannual  variations,  these  biases,  as  well  as  the  mean 
annual  cycle  of  ice  cover,  were  removed.  The  ten-year  average 
annual  cycles  of  A,B  and  C  were  calculated  (Figures  4.27  - 
4.30)  and  subtracted  from  each  of  their  respective  data  fields 
to  produce  anomaly   (difference  from  annual  cycle)  fields  of 
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Figure  4.15  Simulated  thickness  contours  using  prognostic 
ocean  (B)  for  month  76.  Contours  are  cm  of  ice 
thickness . 


89 


T  I  I  I  I  I  r 


1  I  I  1  I  I  I  I  I  I  I  I  I  I  1  I  I  I  I  1  I 


Figure  4.16  Simulated  thickness  contours  using  10-year  mean 
ocean  (C)  for  month  76.  Contours  are  cm  of  ice 
thickness . 
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Figure  4.17  Simulated  thickness  contours  using  prognostic 
ocean  (B)  for  month  93.  Contours  are  cm  of  ice 
thickness . 
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Figure  4.18  Simulated  thickness  contours  using  10-year  mean 
ocean  (C)  for  month  93.  Contours  are  cm  of  ice 
thickness . 
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Figure   4.19 


Time  series  of  total  ice  area  in  region  1. 
Observed  (A)  is  the  solid  line,  prognostic 
ocean  model  (B)  is  the  dashed  line  and  10- 
year  mean  ocean  model  (C)  is  the  dash-dot 
line.  X  axis  is  months  from  Jan  1971  to  Dec 
1980.   Y  axis  is  lO's  of  km^. 
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Figure      4.2  0 


Time  series  of  total  ice  area  in  region  2. 
Observed  (A)  is  the  solid  line,  prognostic 
ocean  model  (B)  is  the  dashed  line  and  10- 
year  mean  ocean  model  (C)  is  the  dash-dot 
line.   X  axis  is  months  from  Jan  1971  to  Dec 


1980 


Y  axis  is  lO's  of  km  . 
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Figure   4.21 


Time  series  of  total  ice  area  in  region  3. 
Observed  (A)  is  the  solid  line,  prognostic 
ocean  model  (B)  is  the  dashed  line  and  10- 
year  mean  ocean  model  (C)  is  the  dash-dot 
line.   X  axis  is  months  from  Jan  1971  to  Dec 


1980 


Y  axis  is  lO's  of  km' 
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Figure   4.2  2 


Time  series  of  total  ice  area  in  region  4. 
Observed  (A)  is  the  solid  line,  prognostic 
ocean  model  (B)  is  the  dashed  line  and  10- 
year  mean  ocean  model  (C)  is  the  dash-dot 
line.   X  axis  is  months  from  Jan  1971  to  Dec 


1980 


Y  axis  is  lO's  of  km. 
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Figure  4.23  Time  series  of  simulated  minus  observed  total  area 


m  region  1. 
dashed  line. 


Dec  1980 


B-A  is  the  solid  line,  C-A  is  the 
X  axis  is  months  from  Jan  1971  to 


Y  axis  is  lO's  of  km' 
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Figure  4.24  Time  series  of  simulated  minus  observed  total  area 

the  solid  line,  C-A  is  the 
is  months  from  Jan  1971  to 
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10 's  of  km' 
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Figure  4.25  Time  series  of  simulated  minus  observed  total  area 
in  region  3.  B-A  is  the  solid  line,  C-A  is  the 
dashed  line.  X  axis  is  months  from  Jan  1971  to 
Dec  1980.   Y  axis  is  lO's  of  km^. 
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Figure  4.26  Time  series  of  simulated  minus  observed  total  area 


m  region  4 
dashed  line 


B-A  is  the  solid  line,  C-A  is  the 


X  axis  IS  months  from  Jan  1971  to 
Dec  1980.   Y  axis  is  10 's  of  km^. 
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Figure  4.27  10-year  mean  annual  cycle  of  total  ice  area  in 
region  1.  Observed  (A)  is  the  solid  line, 
simulated  (B)  is  the  dashed  line  and  simulated  (C) 
is  the  dash-dot  line.   X  axis  is  months  and  Y  axis 

is  lO's  of  km^. 
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Figure  4.28  10-year  mean  annual  cycle  of  total  ice  area  in 
region  2.  Observed  (A)  is  the  solid  line, 
simulated  (B)  is  the  dashed  line  and  simulated  (C) 
is  the  dash-dot  line.  X  axis  is  months  and  Y  axis 
is  10 's  of  km^. 
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Figure  4.29  10-year  mean  annual  cycle  of  total  ice  area  in 
region  3.  Observed  (A)  is  the  solid  line, 
simulated  (B)  is  the  dashed  line  and  simulated  (C) 
is  the  dash-dot  line.  X  axis  is  months  and  Y  axis 
is  IC's  of  km\ 
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Figure  4.30  10-year  mean  annual  cycle  of  total  ice  area  in 
region  4.    Observed   (A)   is  the  solid   line, 
simulated  (B)  is  the  dashed  line  and  simulated  (C) 
is  the  dash-dot  line.   X  axis  is  months  and  Y  axis 

is  10 's  of  kr,  . 
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TABLE  3 

DIFFERENCE  BETWEEN  THE  ABSOLUTE  SUM  OF  SIMULATED 
AND  OBSERVED  ICE  AREA 


REGION 

B  -  A 

C  -  A 

CHANGE 

%  IMPROVEMENT 

1 
2 
3 
4 

.4391 
.1645 
.2945 
.2043 

.4472 
.1729 
.3131 
.2935 

.0081 
.0084 
.0186 
.0892 

2% 

5% 

6% 

43% 

area  in  km^  x  1  x  10® 

ice  concentration  for  A,B  and  C  in  each  region.  The  B  and  C 
anomaly  time  series  are  shown  with  the  A  anomaly  time  series 
for  each  region  in  Figures  4.31  -  4.38.  The  ability  of  the 
B  and  C  anomaly  series  to  follow  the  A  series  was  a  measure 
of  their  ability  to  model  interannual  variability.  A 
quantitative  appreciation  of  this  ability  was  made  in  a 
similar  manner  as  with  the  annual  cycle.  The  A  anomaly  field 
was  subtracted  from  both  B  and  C  and  the  absolute  sum  of  each 
of  the  anomaly  difference  fields  was  then  calculated.  The 
results  are  shown  in  Table  4.  Correlations  were  also 
calculated  between  the  observed  and  modelled  ice  area  fields 
for  both  the  annual  cycle  and  anomaly  time  series.  The 
calculated  correlation  coefficients  are  shown  in  Table  5. 

Several  hundred  contour  plots  and  graphs  were  examined  and 
compared  including  observed  and  simulated  ice  concentration 
and  ice  thickness,  simulated  ocean  heat  flux,  simulated  layer 
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Figure  4.31  Difference  from  mean  annual  cycle  or  anomaly  time 
series  of  total  ice  area  for  region  1.  Observed 
anomaly  is  the  solid  line,  simulated  (B)  anomaly 
is  the  dashed  line.   X  axis  is  months,  Y  axis  is 


10 's  of  km' 
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Figure  4.32  Difference  from  mean  annual  cycle  or  anomaly  time 
series  of  total  ice  area  for  region  1.  Observed 
anomaly  is  the  solid  line,  simulated  (C)  anomaly 
is  the  dashed  line.  X  axis  is  months,  Y  axis  is 
10 's  of  km^ 
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Figure    4.3  3 


Difference  from  mean  annual  cycle  or  anomaly  time 
series  of  total  ice  area  for  region  2 .  Observed 
anom.aly  is  the  solid  line,  simulated  (B)  anomaly 
is  the  dashed  line,   X  axis  is  months,  Y  axis  is 

10 's  of  km^ 
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Figure  4.34  Difference  from  mean  annual  cycle  or  anomaly  time 
series  of  total  ice  area  for  region  2.  Observed 
anomaly  is  the  solid  line,  simulated  (C)  anomaly 
is  the  dashed  line.   X  axis  is  months,  Y  axis  is 

lO's  of  km^ 
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Figure  4.35  Difference  from  mean  annual  cycle  or  anomaly  time 
series  of  total  ice  area  for  region  3.  Observed 
anomaly  is  the  solid  line,  simulated  (B)  anomaly 
is  the  dashed  line.   X  axis  is  months,  Y  axis  is 

10 's  of  km^. 
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Figure  4.36  Difference  from  mean  annual  cycle  or  anomaly  time 
series  of  total  ice  area  for  region  3.   Observed 

anomaly  is  the  solid  line,  simulated  (C)  anomaly 
is  the  dashed  line.  X  axis  is  months,  Y  axis  is 
lO's  of  km^ 
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Figure  4.37  Difference  from  inean  annual  cycle  or  anomaly  time 
series  of  total  ice  area  for  region  4.  Observed 
anomaly  is  the  solid  line,  simulated  (B)  anomaly 
is  the  dashed  line.   X  axis  is  months,  Y  axis  is 

lO's  of  km^. 
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Figure  4.38  Difference  from  mean  annual  cycle  or  anomaly  time 
series  of  total  ice  area  for  region  4 .  Observed 
anomaly  is  the  solid  line,  simulated  (C)  anomaly 
is  the  dashed  line.  X  axis  is  months,  Y  axis  is 
lO's  of  km^ 
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TABLE  4 

DIFFERENCE  BETWEEN  ABSOLUTE  SUM  OF  SIMULATED  AND  OBSERVED 

ICE  AREA  ANOMALY 
(MEAN  ANNUAL  CYCLE  REMOVED) 


REGION 

B  -  A 

C  -  A 

CHANGE 

%  IMPROVEMENT 

1 
2 
3 

4 

.  1476 
.0864 
.  1355 
.0888 

.1818 
.0985 
.1689 
.  1157 

.0342 
.0121 
.0334 
.0268 

23% 
14% 
25% 
30% 

area  in  km  x  1  x  10 


TABLE  5 

CORRELATIONS  OF  ICE  AREA  BETWEEN  SIMULATION  B  AND  OBSERVED 
(CORR  B)  AND  SIMULATION  C  AND  OBSERVED  (CORR  C) 


ICE  AREA  TIME  SERIES 


REG 

STDEV  OBS 

STDEV  B 

STDEV  C 

CORR  B 

CORR  C 

1 
2 
3 
4 

.  3169 
.1434 
.4283 
.2292 

.7776 
.3729 
.6755 
.3667 

.8246 
.3996 
.6910 
.4220 

.9126 
.9271 
.9461 
.9186 

.8852 
.8807 
.9106 
.8332 

ANOMALY  TIKE  SERIES 


REG 

STDEV  OBS 

STDEV  B 

STDEV  C 

CORR  B 

CORR  C 

1 
2 
3 

4 

.1146 
.0590 
.1376 
.0994 

.2220 
.1660 
.1518 
.1311 

.2701 
.1821 
.1677 
.1418 

.4936 
.7761 
.5375 
.6172 

.3963 
.5981 
.3039 
.3332 

Standard  deviation  values  x  1  x  10 

B  simulation  was  with  fully  interactive  ocean  model 
C  simulation  was  with  10-year  mean  cycle  ocean 
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2  velocities  and  the  time  series'  described  above.  On 
conclusion  of  this  review,  the  following  general  observations 
were  made: 

-  The  region  of  transition  from  10/10  to  0/10  ice 
concentration  was  narrower  in  both  B  and  C  than  in  A, 
i.e.,  the  ice  edge  boundary  was  sharper,  more  well  defined 
in  B  and  C  than  in  A. 

-  Both  B  and  C  overpredicted  the  ice  coverage  in  winter, 
especially  in  the  eastern  Barents  Sea  (contained  in  region 
4)  .  B  was  usually  closer  to  the  observed  coverage  than 
C. 

-  Both  B  and  C  underpredicted  the  coverage  in  summer, 
especially  in  the  Kara  Sea  to  the  east  and  the  Beaufort, 
Chukchi  and  East  Siberian  Seas  to  the  west.  B  was  usually 
closer  to  the  observed  coverage  than  C  in  the  eastern 
Arctic  but  the  improvement  in  the  central  and  western 
Arctic  was  minimal. 

-  The  general  shape  of  the  ice  edge  was  reasonably  well 
simulated  by  both  B  and  C.  The  Kara  and  Barents  Seas  were 
notable  exceptions. 

-  The  simulated  surface  (level  2)  currents  matched  the 
observed  annual  mean  currents  (Figure  2.4)  quite  well. 
All  the  major  currents  were  represented  and  their 
velocities  appeared  reasonable;  however,  some  errors  were 
evident  in  the  secondary  currents.  The  West  Spitsbergen 
Current  was  not  represented  and  the  Irminger  Current  was 
flowing  opposite  to  the  observed  mean  direction.  There 
were  significant  differences  between  the  monthly  varying 
fields  and  the  10-year  average  cycle  fields  of  level  2 
ocean  velocities.  The  fields  contained  the  same  current 
features  however  these  features  differed  in  both  strength 
and  position. 

-  The  ocean  heat  flux  into  the  mixed  layer  was  highly 
variable  between  regions  and  changed  dramatically  with 
the  seasons.  Winter  values  were  highest  in  the  Norwegian 
and  Greenland  Seas  where  it  averaged  approximately  180 
W/m  .  Exceptional  cases  of  heat  flux  in  excess  of  375  W/m^ 
were  observed  occasionally  in  a  localized  area  just  south 
of  Svalbard.  Spring  and  summer  heat  fluxes  in  the 
Norwegian  and  Greenland  Seas  were  generally  weakly 
positive  or  negative.  Isolated  summer  heat  flux  values 
less  than  -180  W/m  were  observed  just  north  of  Iceland  in 
some  years.  A  similar  positive  to  negative  ocean  heat 
flux  cycle  was  observed  in  the  western  arctic  although  it 
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was  generally  confined  to  within  400  km  of  the  Alaskan 
shoreline.  The  simulated  ocean  heat  flux  under  the 
central  pack  in  the  Canadian  and  Eurasian  Basins  was 
negligible  throughout  the  year.  Similar  to  the 
velocities,  there  were  significant  differences  between 
the  monthly  varying  fields  and  the  10-year  average  cycle 
fields  in  both  the  strength  and  position  of  the  main 
features. 


-  The  simulated  pack  ice  thicknesses  in  the  central  Arctic 
Basin  were  1-2  m  or  50-100%  less  than  the  average  observed 
values,  particularly  in  summer.  B  thicknesses  tended  to 
be  larger  than  those  for  C.  The  simulated  thickness  data 
indicated  that  some  ridging  was  occurring  in  the  models 
poleward  of  Greenland  and  Ellesmere  Island  as  expected, 
however  the  thickness  values  in  this  region  were  too  low 
by  a  factor  of  three  or  more.  The  observed  thickness 
fields  used  for  comparison  are  shown  in  Figures  4.39  - 
4.42. 

-  Region  1  and  2  had  maximum  ice  areas  limited  by  their 
boundaries.  Region  4's  minimum  ice  area  was  limited  by 
zero  (it  becomes  totally  ice  free).  Region  3  came  close 
to  being  ice  free  especially  from  1978-1980.  These 
minimum  and  maximum  limits  reduced  the  variability  of  the 
total  ice  area;  hov;ever  considerable  yearly  variability 
was  still  evident.  Annual  differences  in  ice  areas  in  the 
observed  data  varied  from  approximately  15%  of  the  mean 
in  region  1  and  2,  to  40%  in  region  3,  to  100%  in  region 
4.  The  simulated  data  had  greater  variance  with 
interannual  ice  area  differences  ranging  from 
approximately  40%  of  the  mean  in  regions  1  and  3,  to  60% 
in  region  2,  to  100%  in  region  4. 

-  The  mean  annual  cycles  of  ice  area  for  A,B  and  C  showed 
that  the  melt/freeze  cycle  in  the  modelled  cases  was 
approximately  in  phase  with  the  observed  annual  cycle  but 
of  higher  amplitude.  C  generally  displayed  the  largest 
amplitude. 

-  The  time  series*  of  ice  area  anomalies  (difference  from 
annual  cycle)  for  both  observed  and  simulated  data  varied 
significantly  from  region  to  region.  Some  similarities 
could  be  found  in  adjacent  regions.  However,  if  they  were 
not  adjacent,  the  correlations  appeared  to  be  very  low. 

-  The  degree  of  monthly  variability  or  "noisiness"  in  the 
anomaly  series  increased  markedly  in  regions  3  and  4 
compared  to  regions  1  and  2. 
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Figure  4.39  Mean  ice  draft  (m)  in  winter  derived  from 
submarine  upward  looking  echo  sounder  data.  Draft 
is  averaged  over  three  months,  centered  on  the 
first  of  r-larch  (from  Bourke  and  Garret,  1987). 
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Figure  4.40  Mean  ice  draft  (in)  in  spring  derived  from 
submarine  upward  looking  echo  sounder  data.  Draft 
is  averaged  over  three  months,  centered  on  the 
first  of  June  (from  Bourke  and  Garret,  1987). 
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Figure  4.41  Mean  ice  draft  (m)  in  summer  derived  from 
submarine  upv;ard  looking  echo  sounder  data.  Draft 
is  averaged  over  three  months,  centered  on  the 
first  of  September  (from  Bourke  and  Garret,  1987)  . 
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Figure  4.42  Mean  ice  draft  (m)  in  fall  derived  from  submarine 
upward  looking  echo  sounder  data.  Draft  is 
averaged  over  three  months,  centered  on  the  first 
of  December  (from  Bourke  and  Garret,  1987). 
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Figure  4.43  10-year  mean  ocean  heat  flux  through  the  30 
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The  B  anomaly  time  series  was  a  better  match  to  the 
observed  (A)  anomaly  time  series  than  C.  There  were  time 
periods  when  both  B  and  C  were  significantly  different 
from  A  but  with  the  same  sign.  There  are  also  times  when 
they  both  differed  from  the  A  anomalies  but  with  opposite 
signs.  No  clear  pattern  was  evident  of  a  bias  to  either 
sign. 

Considering  the  entire  12  0  month  time  period,  model  runs 
using  the  interactive  ocean  (B)  matched  the  observed  data 
better  than  the  average  ocean  case  (C)  in  all  regions. 
This  was  true  for  both  total  ice  area  and  interannual 
variability.  Both  the  absolute  sums  of  total  differential 
area  and  the  correlation  coefficients  for  each  region  were 
better  for  B  than  for  C,  particularly  in  region  4. 
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C.   DISCUSSION 

The  excessive  melting  which  occured  in  regions  1  and  2  was 
a  common  feature  of  both  the  prescribed  ocean  model  and  the 
prognostic  ocean  model.  This  suggests  that  an  error  or 
weakness  common  to  both  model  formulations  may  be  responsible. 
One  possible  weakness  is  the  quite  simple  representation  of 
the  mixing  processes  within  the  surface  layer  and  the  rather 
large  vertical  resolution  of  the  upper  500  m.  Contours  of 
ocean  heat  flux  (e.g.,  Figure  4.43)  show  large  positive  values 
in  summer  when  the  melting  rate  is  highest.  The  large  amount 
of  melting  in  this  season  would  be  expected  to  produce  a 
surface  cap  of  fresh  water.  The  resultant  stratification 
would  insulate  the  mixed  layer  from  the  heat  below  and  reduce 
this  flux.  It  appears  that  this  mechanism  is  not  being 
represented  properly  as  excessive  heat  continues  to  enter  the 
mixed  layer  until  late  summer.  Alternatively,  the  large 
lateral  diffusion  coefficients  which  are  necessary  for 
numerical  stability  may  be  mixing  far  too  much  heat  throughout 
the  intermediate  layer.  In  this  case,  even  if  the  vertical 
advection  is  reduced  by  stratification,  there  may  still  be 
enough  heat  available  to  cause  the  excessive  summer  ice 
retreat. 

A  second  possibility  becomes  evident  after  sequentially 
viewing  the  ocean  heat  flux  contours  over  several  months.  A 
close  examination  of  those  areas  which  appear  most  in  error 
indicates  several  "hot"  spots  in  winter  and  spring  which 
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correspond  to  some  of  the  locations  of  prescribed  inflow. 
These  same  "hot"  spots  become  "cold"  spots  during  summer  and 
fall.  The  constant  prescribed  inflow  temperatures  at  the 
Bering  Strait  and  Faeroe-Shetland  Channel  contrast  with  the 
simulated  ocean  temperatures  produced  by  the  model.  The 
simulated  temperatures  vary  monthly  following  the  annual  cycle 
of  SST  resulting  in  a  cycle  of  anomalous  temperatures  at  the 
inflow  positions.  Although  the  temperature  of  the  river 
inflow  is  set  to  match  the  ocean  temperature,  similar  "hot" 
spots  occur  at  the  river  mouths  too,  particularly  for  the 
Mackenzie,  Kolyma  and  Lena  Rivers.  The  most  probable  causes 
for  this  appear  to  be  geostrophic  adjustment  to  the  negative 
salt  fluxes  and/or  bottom  topography  interaction.  The  heat 
from  these  hot  spots  is  advected  by  the  currents,  spreading 
their  effect  and  increasing  the  spring  melting  over  a 
considerable  area.  The  western  boundary  areas  appear  to  be 
a  good  example  of  this. 

The  simulated  v;inter  ice  cover,  although  not  an  exact 
match  to  the  observed,  is  still  a  closer  match  than  the  summer 
condition.  This  would  seem  reasonable  because  the  large 
freezing  rates  and  surface  cooling  in  the  winter  induce  deep 
vertical  mixing.  This  process  can  be  better  represented  by 
the  simple  density  mixing  scheme  used  in  this  model  than  the 
shallow  stratification  process  which  occurs  in  summer. 
Realistic  amounts  of  heat  are  brought  to  the  surface  and  the 
simulated  ice  edge  matches  the  observed  quite  well.    The 
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improvement  provided  by  the  prognostic  ocean  in  regions  3  and 
4  is  consistent  because,  if  we  consider  that  the  winter 
vertical  convection  is  reasonably  well  accounted  for,  an 
improvement  in  the  interannually  varying  currents  and  ocean 
heat  flux  should  have  a  direct  improvement  on  the  ice  cover 
simulation. 

D.   CONCLUSIONS 

Inclusion  of  the  interactive  ocean  produced  an  improved 
accuracy  of  simulation  of  the  annual  cycle  of  total  ice  area, 
particularly  in  region  4 .  The  model  which  included  the 
interactive,  prognostic  ocean  represented  the  variable  ocean 
forcing  (heat  and  momentum)  much  better  than  the  average  ocean 
model.  This  suggests  that  variable  ocean  forcing  is  important 
to  the  evolution  of  the  ice  cover  in  the  eastern  Arctic.  By 
contrast,  in  regions  1  and  2,  the  interactive  model  did  not 
improve  the  representation  of  the  total  ice  area 
significantly.  In  these  regions,  it  may  be  the  atmospheric 
forcing,  the  ice  rheology  or  as  discussed  above,  mechanisms 
controlling  mixing  of  heat  in  the  upper  ocean  layers,  that 
dominate  the  evolution  of  total  ice  area.  These  are 
parameters  which  were  not  changed  between  experiments. 
Alternatively,  there  could  be  errors  in  the  ocean  forcing 
which  produces  a  bias  in  the  ice  field  in  both  the  variable 
ocean  and  the  average  ocean. 
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Inclusion  of  the  interactive  model  produced  significant 
improvements  in  simulating  the  interannually  varying  ice  field 
in  all  regions  of  the  Arctic.  This  indicates  that  the 
interannual  changes  in  the  ocean  forcing  are  a  major  influence 
on  the  interannual  variability  of  the  ice  field  over  the 
entire  Arctic  basin.  The  importance  of  including  an 
interactive  ocean  in  seasonal  ice  prediction  models  is 
therefore  demonstrated. 

Simulated  ice  thickness  has  not  been  improved 
substantially  based  on  the  data  presented  by  Bourke  and  Garret 
(1987).  Again,  as  for  total  ice  area  in  the  western  Arctic, 
modifications  to  the  way  atr.ospheric  forcing  is  handled  (e.g., 
albedo,  surface  drag)  or  to  the  ice  rheology  may  be  required 
to  improve  this  feature  of  the  ice  cover. 
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V.   ICE  STRENGTH  EXPERIMENT 

One  of  the  most  obvious  inaccuracies  in  the  simulated  ice 
cover  produced  by  the  linked  model  was  the  ice  thickness  as 
discussed  in  Chapter  IV.  In  the  polar  pack  region,  simulated 
thicknesses  were  approximately  one  half  the  observed  values. 
Thicknesses  in  the  ice  convergence  regions  north  of  Greenland 
and  the  Canadian  archipelago  were  too  small  by  a  factor  of 
three  or  more. 

Heavy  ice  buildup  in  the  convergence  regions  is  the  result 
of  the  large  amount  of  ridging  and  rafting  which  occurs  there. 
The  poor  representation  of  this  feature  indicated  that  the  ice 
rheology  used  was  not  portraying  these  mechanisms  properly. 
Hibler  (1979,  1980)  was  able  to  simulate  reasonable  ice 
thicknesses  in  the  strong  convergence  areas  along  the  north 
coast  of  Ellesmere  Island,  but  noted  that  the  simulated  ice 
buildup  in  this  region  and  overall  spatial  ice  thickness 
variations  were  dependent  on  the  ice  strength  value  used.  In 
the  follow-on  work  by  Hibler  and  Bryan  (1987) ,  a  greater  ice 
strength  was  used  to  provide  more  accurate  ice  velocities  at 
the  larger  grid  scale.  However,  the  ice  thickness  in  the 
convergence  regions  was  considerably  reduced,  appearing  quite 
similar  to  the  initial  results  found  in  this  work. 

Examination  of  the  ice  strength  value  used  in  the  model 
here  indicated  that  the  ice  strength  parameter,   P  ,  was 
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relatively  large.  This  made  the  ice  pack  quite  stiff, 
reducing  the  amount  of  ridging  in  areas  of  ice  convergence. 
The  ice  rheology  used  in  this  model  was  the  same  as  in  the 
Semtner  (1987)  model.  However,  it  should  be  noted  that  the 
strength  parameter  actually  used  by  Semtner  (1987)  differed 
from  the  description  contained  in  the  Appendix  to  that 
reference.   The  P  actually  used  was: 

P*  =  h  X  3  X  10^  (N/m^) 

where  h  is  the  ice  thickness.  This  was  30  times  larger  than 
the  value  noted  in  the  appendix.  Semtner  indicated  (private 
communication)  that  this  larger  value  was  required  to  prevent 
the  numerical  instability  which  had  occurred  in  earlier  runs 
using  a  smaller  P  value. 

This  experiment  investigated  the  effects  of  reducing  P  . 
The  objective  was  to  determine  the  best  value  for  producing 
an  accurate  simulation  of  the  ice  thickness  and  thickness 
distribution. 

A.   EXPERIMENTAL  SETUP 

The  fully  prognostic  ice-ocean  model  (B)  was  used  for  this 
experiment.  Initial  runs  with  varying  P  values  indicated 
that  changes  to  the  simulated  ice  fields  were  dramatic. 
Indeed,  changes  were  evident  after  only  a  single  month  of 
integration.  Therefore,  single  year  integrations  were 
considered  sufficient  for  the  sensitivity  study  vice  a  ten 
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year  integration  period.   The  year  1974  was  selected  for  the 
runs  as  it  was  a  fairly  average  ice  cover  year. 

Seven  runs  were  conducted.  P  values  were  set  as  follows: 


-  P*  =  h  X  3  X  10^ 

-  P*  =  h  X  1  X  10^ 


-  P*  =  h  X  5  X  10^ 


-  P*  =  h  X  3  X  lo"^ 


-  P*  =  h  X  1  X  10^ 


-  P*  =  h  X  5  X  10^ 


-  P*  =  h  X  1  X  10^ 


Once  an  optimum  P  was  determined,  the  full  ten  year 
integration  was  conducted  to  produce  ice  concentration  fields 
for  comparison  with  the  observed  data. 

B.   RESULTS 

Numerical  instability  occurred  in  this  model  at  the 
reduced  P  values  similar  to  Semtner's  work  (private 
coFimunication)  .  The  instability  originated  at  a  grid  point 
in  the  EGC,  adjacent  to  where  the  grid  representation  of  the 
Greenland  coastline  was  very  jagged.  The  jagged  coastline  was 
a  result  of  the  model's  limited  resolution.  Ocean  current 
velocities  grew  exponentially,  producing  artificially  high  ice 
velocities  and  excessive  kinetic  energy.  The  temperature  and 
salinity  values  also  became  unreasonable.  Reduction  of  the 
model's  numerical  time  steps  by  a  factor  of  two  did  not  solve 
the  problem.   Numerical  stability  was  eventually  obtained  by 


129 


imposing  velocity,  temperature  and  salinity  limits.    The 
limitations  on  these  prognostic  variables  were  as  follows; 
-75.0  <  UICE,VICE  <  +75.0    (cm/sec) 
-2.0  <  TSURF  <  15.0    (deg  centigrade) 
0  <  SALINITY  <  50    (ppt) . 

Figures  5.1  -  5.6  show  the  simulated  ice  thickness 
contours  for  August  1974,  using  the  P  values  from  runs  2 
through  7,  respectively.  In  general,  as  P  was  reduced,  the 
entire  ice  pack  tended  to  shift  toward  the  Fram  Strait  and 
thicknesses  within  the  pack  increased.  The  convergence  along 
the  north  coast  of  Greenland  became  more  intense  eventually 
producing  a  reasonable  simulation  of  the  ridging  and  rafting 
processes.  However  the  simulated  region  of  maximum 
thicknesses  was  still  thinner  and  was  located  further  poleward 
than  the  average  thickness  fields  in  Bourke  and  Garrett  (1987) 
(See  Figures  4.39  -  4.42).  The  tendency  for  the  model  to 
produce  too  much  meltback  of  the  ice  edge  in  summer  was 
exacerbated,  even  with  the  increased  ice  thicknesses.  This 
was  especially  evident  in  the  Bering  and  East  Siberian  Seas. 

The  thickness  gradients  became  increasingly  stronger, 
irregular  and  spotty  as  the  ice  strength  was  set  to  the  lower 
values.  The  two  lowest  P  values  produced  thickness  contours 
which  appeared  too  "noisy"  to  be  realistic. 

A  qualitative  assessment  of  all  the  thickness  contour 
plots  and  comparison  with  the  Bourke  and  Garrett  (1987)  data 
indicated  that  P*  =  h  x  1  x  10^  provided  the  most  reasonable 
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Figure  5.1   Simulated  ice  thickness  contours  (cm)  for  August, 
1974.     P*  =  h  X  1  X  10^ 
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Figure  5.2   Simulated  ice  thickness  contours  (cm)  for  August, 
1974  .     P*  =  h  X  5  X  10\ 
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Figure  5.3   Simulated  ice  thickness  contours  (cm)  for  August, 
1974  .     P*  =  h  X  3  X  10\ 
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Figure  5.4   Simulated  ice  thickness  contours  (cm)  for  August, 
1974  .     P*  =  h  X  1  X  10\ 
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Figure  5.5   Simulated  ice  thickness  contours  (cm)  for  August, 
1974.     P*  =  h  X  5  X  10^. 
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Figure  5.6   Simulated  ice  thickness  contours  (cm)  for  August, 
1974  .     p'  =  h  X  1  X  10^. 
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simulated  overall  ice  thickness  and  thickness  distributions. 
This  value  is  somewhat  larger  than  that  used  in  Hibler 
(1979,1980)  (P*  =  f(h)  X  5  X  10^)  and  smaller  than  that  used 
in  Hibler  and  Bryan  (1987)  (P*  =  f(h)  x  2.75  x  10^).  Both  of 
those  were  chosen  by  matching  simulated  drift  rates  with 
observed  data.  The  reason  for  a  different  optimum  value  for 
P  is  probably  due  to  the  higher  spatial  resolution  and 
monthly  averaged  daily  atmospheric  forcing  used  in  this  model 
vice  the  observed  daily  forcing  used  in  the  Hibler  models. 

The  ten-year  integration  was  run  with  the  optimum  P  and 
a  fully  interactive  ocean  to  produce  ten  years  of  simulated 
ice  concentration  fields.  Correlations  with  the  observed  data 
for  both  the  annual  cycle  and  interannual  variations  for  all 
four  regions  were  calculated.  The  correlations  using  the 
reduced  P  model  were  considerably  lower  than  when  the 
experiment  1  model  with  a  large  P  value  was  used  (Table  6)  . 

C.   DISCUSSION 

The  linked  ice-ocean  model  with  a  prognostic  ocean 
displayed  notable  sensitivity  to  changes  in  the  ice  strength 
factor  of  the  ice  rheology.  Reducing  the  P  value  by  a  factor 
of  3  0  improved  the  simulated  ice  thickness  and  distribution 
significantly.  However,  further  reductions  in  ice  strength 
appeared  to  introduce  excessive  "noise"  in  these  fields.  The 
contour  lines  became  very  rough  and  erratic  with  numerous 
"bullseyes"  or  spot  irregularities.  This  was  probably  due  to 
the  diminished  ability  of  the  ice  to  sustain  internal  stress. 
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TABLE  6 

CORRELATIONS  OF  ICE  AREA  BETWEEN  SIMULATION  B  AND  OBSERVED 
(CORR  B)  AND  SIMULATION  D  AND  OBSERVED  (CORR  D) 


AREA  TIME  SERIES 


REG 

STDEV  D 

CORR  B 

CORR  C 

CORR  D 

1 
2 
3 
4 

.8160 
.3882 
.6709 
.4394 

.9126 
.9271 
.9461 
.9186 

.8852 
.8807 
.9106 
.8332 

.8610 
.8866 
.9322 
.8623 

ANOMALY  TIME  SERIIES 


REG 

STDEV  D 

CORR  B 

CORR  C 

CORR  D 

1 
2 
3 
4 

.1697 
.1499 
.1748 
.1788 

.4936 
.7761 
.5375 
.6172 

.3963 

.5981 
.3039 
.3332 

.4294 
.6292 
.4536 
.4354 

Standard  deviation  values  x  1  x  10 

B  simulation  was  with  fully  interactive  ocean  model 
C  simulation  was  with  10-year  mean  cycle  ocean 
D  simulation  was  with  the  reduced  P 
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As  the  ice  strength  was  reduced,  the  ice  approached  a  free 
drift  condition  with  little  resistance  to  either  shearing  or 
compressive  forces.  In  this  condition,  the  ice  could  respond 
rapidly  and  easily  to  any  short  term  or  small  scale  (a  few 
hundred  km)  dynamic  forcing  resulting  in  the  irregular 
appearance  of  the  thickness  contours  at  the  lowest  P  values. 
Since  the  atmospheric  forcing  was  smoothed,  these  shorter  or 
small  scale  forcing  events  were  probably  from  the  ocean. 

An  explanation  for  the  increased  ice  thicknesses  when  a 
moderate  reduction  in  ice  strength  was  introduced  is  that  by 
reducing  the  ice  strength,  the  long  term  dynamic  influences 
of  the  Transpolar  Drift  Stream  and  atmosphere  were  able  to 
push  the  ice  further  and  faster  to  the  southeast.  The 
resultant  ice  divergence  from  the  Bering  and  adjacent  seas 
created  large  expanses  of  open  water.  Large  heat  losses  to 
the  atmosphere  associated  with  this  open  water  and  cold  air 
temperatures  in  fall  and  winter  encouraged  rapid  ice 
production.  The  large  volumes  of  new  ice  were  then  advected 
to  the  southeast  where  at  the  lower  strength  levels  they  were 
able  to  compress  and  undergo  the  appropriate  ridging  and 
rafting  processes.  This  produced  the  large  areas  of  thin  ice 
on  the  Siberian  shelves,  increased  the  average  ice  pack 
thicknesses  and  resulted  in  the  stronger  thickness  gradients 
evident  in  the  simulated  ice  thickness  fields. 

The  ice  thickness  and  distribution  fields  from  the  reduced 

* 
P   model   appeared   to   be   so   much   improved,   that   the 
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correlations  of  ice  concentration  with  the  observed  fields 
were  expected  to  improve  also.  The  subsequent  determination 
that  these  correlations  were  actually  degraded  was  surprising. 
Ice  thickness  and  distribution  appeared  to  improve 
substantially  within  the  ice  pack.  However  the  amount  of 
ocean  area  which  experienced  a  total  loss  of  ice  actually 
increased,  particularly  in  the  Bering  and  East  Siberian  Seas. 
Since  this  exaggerated  the  already  excessive  ice  edge  retreat 
in  summer,  and  no  compensatory  improvement  in  the  ice  edge 
position  was  produced  in  the  other  seasons,  the  correlations 
with  the  observed  data  decreased. 

The  increased  ice  edge  retreat  and  areas  of  open  water  can 
be  explained  as  follows.  The  large  heat  loss  over  the  open 
water  and  the  strong  positive  salt  flux  from  the  large 
freezing  rate  would  promote  increased  convective  overturning 
of  the  ocean  in  the  Siberian  shelf  region.  Strong  convection 
would  break  through  the  halocline  which  serves  to  insulate  the 
surface  waters  from  the  warmer  waters  below.  Oceanic  heat 
flux  to  the  ocean  surface  would  increase,  slowing  the  freezing 
rate  as  winter  progressed  and  keeping  the  ice  cover  relatively 
thin  in  those  regions.  Once  the  increased  solar  flux  in 
spring  produced  a  net  positive  heat  flux  at  the  ocean  surface, 
the  thin  ice  would  rapidly  melt  and  extensive  areas  of  open 
water  would  form. 

Hibler   and   Bryan    (1987)    describe   the   complex 
interactions  and  feedback  mechanisms  possible,  particularly 
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in  the  marginal  ice  zone,  when  the  effects  of  dynamics  and  an 
interactive  ocean  are  included.  They  note  that  a  fine  balance 
of  ice  processes  appears  to  exist.  The  initial  freezing 
causes  intense  overturning.  The  resultant  deeply  mixed 
surface  waters  can  withstand  considerable  amounts  of  melting 
at  a  later  time  without  becoming  stratified.  This  description 
is  consistent  with  the  process  suggested  here. 

D.   CONCLUSION 

An  appropriate  function  to  use  for  ice  strength  in  the  ice 
rheology  of  this  model  appears  to  be: 

P*  =  h  X  1  X  10^ 

The  ice  thickness  and  distribution  fields  produced  using  this 
function  were  considerably  more  realistic  in  comparison  with 
available  observations  than  those  produced  using  a  greater  ice 
strength.  The  simulated  average  pack  thickness  was  doubled 
and  strong  ice  buildup  occurred  in  a  manner  and  in  areas 
consistent  with  observed  data.  A  possible  exception  was  the 
region  immediately  poleward  of  Greenland  and  the  Canadian 
Archipelago.  Ice  thicknesses  there  were  reduced,  apparently 
due  to  advection  because  there  did  not  appear  to  be 
significant  ocean  heat  flux.  Smaller  ice  strength  values 
produced  thickness  fields  which  were  considered  too  noisy. 

Both  the  annual  cycle  and  interannual  variations  of  ice 
concentration  produced  using  the  reduced  ice  strength  were 
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degraded,  primarily  because  the  total  ocean  area  which  became 
ice  free  in  summer  increased  from  the  initial  model  runs.  It 
is  proposed  that  the  reason  for  the  increased  area  of  open 
ocean  in  summer  is  an  increase  in  the  ocean  heat  flux.  This 
occurs  as  a  result  of  the  strong  convective  penetration  of  the 
halocline  resulting  from  greater  heat  losses  to  the  atmosphere 
from  regions  where  the  ice  has  been  thinned  or  removed  and 
increased  salt  fluxes  into  the  ocean  caused  by  the  larger 
freezing  rates  in  the  same  regions. 

These  results  indicate  that  the  simulated  heat  flux 
provided  by  the  ocean  appears  to  be  a  dominant  factor 
controlling  the  differences  between  simulated  and  observed  ice 
edge  position.  Inclusion  of  the  interactive  ocean  in  the 
linked  ice-ocean  model  provided  a  much  more  realistic 
representation  of  this  flux.  Consequently,  simulated  ice 
fields  were  much  improved  over  the  model  which  used  only  an 
average  ocean  condition.  However,  the  inaccuracies  that 
remain  in  the  representation  of  ocean  heat  flux  appear  to  be 
very  important. 
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VI.   ALBEDO  EXPERIMENT 

The  excessive  meltback  of  the  ice  edge  in  summer  was  a 
common  feature  of  previous  model  runs  in  this  work.  Similar 
results  have  been  reported  by  other  ice  modellers  as  well 
(e.g.,  Semtner,  1987;  Hibler  and  Walsh,  1982).  The  ice 
strength  experiment  indicated  that  this  may  have  been  due  to 
excessive  simulated  oceanic  heat  flux  in  the  areas  of  severe 
melting.  Alternatively,  excessive  atmospheric  heat  flux  into 
the  ice  cover  could  also  induce  an  exaggerated  meltback  of  the 
ice  edge.  One  of  the  most  important  elements  determining  the 
atmospheric  heat  flux  into  the  ice  is  surface  albedo. 

Shine  and  Henderson-Sellers  (1985)  noted  the  differences 
in  albedo  representations  used  by  several  authors  (e.g., 
Parkinson  and  Washington,  1979;  Hibler  and  Walsh,  1982;  Manabe 
and  Stouffer,  1980) .  They  showed  that  a  thermodynamic  ice 
model  very  similar  to  the  one  incorporated  into  this  model  was 
quite  sensitive  to  changes  in  the  surface  albedo.  Large 
differences  between  various  model  thickness  predictions  could 
be  explained  by  these  albedo  differences. 

The  Shine  and  Henderson-Sellers  (1985)  study  used  a  model 
which  did  not  include  ice  dynamics  nor  an  interactive  ocean. 
As  noted  in  Chapter  V,  inclusion  of  these  in  a  linked  model 
permits  a  much  more  realistic  representation  of  the  complex 
mixing  and  feedback  processes  which  occur  during  freezing  and 
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melting.  The  response  of  the  more  elaborate  model  used  in 
this  study  to  surface  albedo  changes  was  unknown  and  therefore 
was  considered  worthwhile  for  further  investigation. 

A.   EXPERIMENTAL  SETUP 

The  represention  of  surface  albedo  used  in  this  model  was 
described  in  Chapter  III.  Some  allowance  was  made  for  the 
differences  in  albedo  between  ice  and  snow  and  for  the  changes 
which  occur  when  melting  commences  (see  Chapter  III)  ;  however, 
the  representation  was  very  simple. 

Ross  and  Walsh  (1987)  used  a  more  elaborate  surface  albedo 
representation  in  conjunction  with  a  dynamic-thermodynamic 
model  to  simulate  surface  albedo  in  the  Arctic.  The  albedo 
representation  was  a  linear  function  of  air  temperature  which 
also  accounted  for  the  differences  between  ice  and  snow.  The 
simulated  albedos  produced  with  that  model  matched  the 
satellite-derived  estimates  by  Robinson  et  al.,  (1986)  very 
well.  In  view  of  the  importance  placed  on  an  accurate 
representation  of  the  surface  albedo  by  Shine  and  Henderson- 
Sellers  (1985),  the  albedo  representation  of  Ross  and  Walsh 
(1987)  was  incorporated  into  the  model  used  in  this  work  as 
follows: 

-  a^^ow  =0.80  T3,,  <  -5.0  deg  C 

-  ^snow  =  0.65  -  0.03(T3,^  )    -5.0  <  T^,,  <  0.0  deg  C 

-  ^snow  =  0.65  T,,,  =  0.0  deg  C 

-  a,,^   =0.65  T^^,  <  0.0  deg  C 
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-  a,.^g   =  0.65  -  0.04(T3,-,)        0.0  <  Tg,-,  <  5.0  deg  C    ^ 

-  ^^ce      =  0.45  T3^,  >  5.0  deg  C 

-  a,3,er  =  0.10 

where  "a"  is  the  albedo  value,  T^^^  is  the  surface  temperature 
of  the  ice  or  snow  and  T^,,.  is  the  air  temperature  above  the 
ice  or  snow  surface.  The  fraction  of  solar  radiation  absorbed 
internally  by  the  ice  was  0.35  as  recommended  by  Ross  and 
Walsh  (1987). 

A  ten  year  integration  was  conducted  using  the  new  albedo 
model  and  radiation  absorption  constant.  The  resultant  ice 
concentration  fields  were  again  correlated  with  the  observed 
data  as  in  the  previous  experiments.  The  ten-year 
correlations  of  simulated  annual  ice  cycle  and  interannual 
variations  are  shown  in  Table  7.  The  correlation  coefficients 
improved  nominally  in  regions  2  and  4  and  worsened  nominally 
in  regions  1  and  3. 

Comparison  of  simulated  ice  concentration  contours  and  ice 
thickness  contours  from  model  runs  before  and  after  the  albedo 
changes  showed  only  trivial  differences.  It  appeared  that  the 
full  model  with  ice  dynamics  and  an  interactive  ocean  was 
relatively  insensitive  to  changes  in  albedo.  This  finding  was 
a  clear  contrast  to  previous  work  using  simpler  ice  models 
(e.g.,  Shine  and  Henderson-Sellers,  1985;  Hibler,  1980). 


note  that  this  is  a  correction  to  Ross  and  Walsh  (1987) 
to  make  the  albedo  a  linear  function  of  temperature  as 
described  in  the  text. 
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TABLE  7 

CORRELATIONS  OF  ICE  AREA  BETWEEN  SIMULATION  B  AND  OBSERVED 
(CORR  B)  AND  SIMULATION  D  AND  OBSERVED  (CORR  D) 


ICE  AREA  TIME  SERIES 


REG 

STDEV  E 

CORR  B 

CORR  D 

CORR  E 

1 
2 
3 
4 

.8066 

.3778 
.6647 
.4387 

.9126 
.9271 
.9461 
.9186 

.8610 
.8866 
.9322 
.8623 

.8606 
.8870 
.9321 
.8648 

ANOMALY  TIME  SERIES 


REG 

STDEV  E 

CORR  B 

CORR  D 

CORR  E 

1 
2 

3 

4 

.1715 
.1491 
.1767 
.  1798 

.4936 
.7761 
.5375 
.6172 

.4294 
.6292 
.4536 
.4354 

.4277 
.6328 
.4522 
.4556 

Standard  deviation  values  x  1  x  10  . 

B  simulation  was  with  fully  interactive  ocean  model. 

C  simulation  was  with  10-year  mean  cycle  ocean. 

D  simulation  was  with  the  reduced  P  . 

E  simulation  was  with  the  improved  albedo  representation, 
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Several  sensitivity  runs  were  made  to  further  explore  this 
observation.  The  range  of  the  albedos  was  systematically 
increased  until  values  well  outside  the  normal  range  were 
prescribed.  This  was  done  to  determine  just  how  insensitive 
the  model  was  to  albedo  modifications.  These  runs  used  a 
single  year  of  integration  (1977)  in  the  same  manner  as  for 
the  ice  strength  sensitivity  experiment.  The  surface  albedo 
representations  examined  were  as  follows: 

-  Al.   Albedo  code  as  in  Semtner  (1987). 

-  A2 .   New  albedo  code  as  described  above  and  used  in  Ross 

and  Walsh  (1987)  . 

-  A3.   All  ice  and  snow  albedos  in  the  new  albedo  code 

increased  by  0.05. 

-  A4 .   All  ice  and  snow  albedos  in  the  new  albedo  code 

increased  by  0.1  and  open  water  albedo  increased  to 
0.15. 

-  A5 .   Ice  albedo  held  constant  at  0.65,  snow  albedo  held 

constant  at  0.80  and  open  water  albedo  at  0.15. 

-  A6.   Ice  albedo  held  constant  at  0.85,  snow  albedo  held 

constant  at  0.90  and  open  water  albedo  at  0.15. 

-  A7 .   Ice  albedo  held  constant  at  0.35,  snow  albedo  held 

constant  at  0.4  0  and  open  water  albedo  at  0.10. 

Several  other  runs  were  conducted  using  the  high  and  low 
albedo  representations  of  runs  A6  and  A7  but  with  different 
portions  of  the  thermodynamic  and  dynamic  forcing  in  the 
linked  model  removed.  This  was  done  to  determine  which  of  the 
features  of  the  fully  linked  model  reduced  its  sensitivity  to 
albedo  changes  as  compared  to  previous  thermodynamic  ice 
models.  This  work  was  done  in  conjunction  with  the  objectives 


147 


of  Chapter  VII.  For  reference,  a  brief  description  of  the 
various  dynamic  model  variations  are  noted  below.  Further 
descriptions  of  the  various  runs  Dl  thru  D8  and  the  rationale 
for  their  selection  are  contained  in  Chapter  VII. 

-  Dl.   Base  model  using  the  10-year  average  ocean.   Start 

from  Dec  1976  and  use  1977  forcing  over  a  single 
integration  year. 

-  D2 .   Constant  ocean  heat  flux  into  mixed  layer  (2.0  x  10' 

^  deg  C/sec,  equivalent  to  approximately  25  W  m'^)  . 

-  D3 .   No  ocean  heat  flux. 

-  D4 .   No  ocean  current  stress  (level  two  (40  m)  ocean 

velocities  =  0.0  m/s) . 

-  D5 .   Ocean  heat  flux  =  2.0  x  10'  and  no  ocean  current 

stress. 

-  06.   No  wind  stress  (surface  drag  coefficient  =  0.0). 

-  D7 .   No  wind  or  ocean  current  stress. 

-  D8 .   No  wind  stress,  no  ocean  current  stress  and  no  ocean 

heat  flux. 

In  cases  02  through  08,  the  description  for  each  run  indicates 
the  changes  from  the  Base  Model  (run  01) . 

B.   RESULTS 

The  ice  thickness  distributions  for  September  1977  from 
runs  A6  and  A7  are  shown  in  Figures  6.1  and  6.2.  These  two 
runs  had  the  largest  differences  in  albedos  and  September  was 
the  month  which  showed  the  greatest  differences  in  ice 
thicknesses.  The  differences  between  runs  A6  and  A7  (and  runs 
A1-A5  as  well)  were  small.  The  ice  cover  did  in  general 
increase  when   the  albedo  was   changed  from  a   low  to  a  high 
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Fiqure  6.1   Simulated   ice   thickness   contours   (cm)   for 


September,  1977.   High  albedo  case  (run  AS; 
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Figure  6.2   Simulated   ice   thickness   contours   (cm] 
September,  1977.   Low  albedo  case  (run  A7) 
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value;  however,  the  ice  edge  advance  rarely  exceeded  50  km 
seaward.  In  fact,  there  were  numerous  examples  in  which  the 
position  of  the  ice  edge  and  contours  within  the  high 
thickness  gradient  regions  were  indistinguishable  between  the 
high  and  low  albedo  cases  except  for  the  smallest  of 
irregularities . 

Differences  between  simulations  using  modified  albedos 
became  more  noticeable  in  runs  D1-D8.  The  common  differences 
were  a  southward  extension  of  the  marginal  ice  zone  (MIZ)  all 
around  the  Arctic  Basin,  a  southward  extension  of  the  EGC, 
increased  ice  in  the  Kara  Sea  and  in  several  cases  an  increase 
in  ice  poleward  of  Spitsbergen.  Summer  central  ice  pack 
thicknesses  decreased  very  slightly  (approximately  2%)  in 
several  of  the  high  albedo  runs  presumably  because  the  thicker 
ice  around  the  margins  reduced  the  amount  of  ice  comprer,sion 
in  the  central  pack.  Differences  between  high  and  low  albedo 
simulations  were  largest  in  the  spring  during  periods  of  large 
scale  melting.  Once  the  melt  was  in  full  progress  (and  during 
the  subsequent  freezing  period) ,  the  ice  fields  evolved  almost 
identically. 

Figures  6.3  -  6.8  show  the  ice  thickness  distributions  in 
July  1977  for  the  high  and  low  albedo  conditions  from  runs  D6, 
D7  and  D8 .  These  were  the  runs  in  which  wind  stress  was 
reduced  to  zero.  In  these  cases  much  more  ice  was  evident  in 
the  eastern  Arctic  during  summer  and  fall  using  the  large 
albedo  values.   Increased  southward  extension  of  the  MIZ  was 
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Figure  6.3   Simulated  ice  thickness  contours  (cm)  for  July, 
1977.   High  albedo  case  (run  D6) . 
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Figure  6.4   Simulated  ice  thickness  contours  (cm)  for  July, 
1977.   Low  albedo  case  (run  D6) . 
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Figure  6.5   Simulated  ice  thickness  contours  (cm)  for  July, 
1977.   High  albedo  case  (run  D7) . 
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Figure  6.6   Siinulated  ice  thickness  contours  (cm)  for  July, 
1977,   Low  albedo  case  (run  D7). 
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Figure  6.7   Simulated  ice  thickness  contours  (cm)  for  July, 
1977.   High  albedo  case  (run  D8). 
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Figure  6.8   Simulated  ice  thickness  contours  (cm)  for  July, 
1977.   Low  albedo  case  (run  D8). 
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also  observed  in  the  western  Arctic  but  relative  to  the 
eastern  Arctic  increase,  this  change  was  relatively  small. 

In  contrast  to  the  small  differences  in  ice  cover  between 
simulations  when  the  albedo  was  changed,  the  differences 
between  simulations  when  the  dynamic  or  thermodynamic  forcing 
was  changed  were  dramatic.  These  effects  will  be  examined  in 
Chapter  VII. 

C.   DISCUSSION 

Prior  to  this  experiment,  it  was  tempting  to  assume  that 
the  simple  albedo  representation  used  by  Semtner  (1987)  and 
initially  used  here  provided  surface  albedo  values  which  were 
too  low.  For  a  model  sensitive  to  albedo  changes,  this  would 
cause  earlier  and  increased  melting  and  perhaps  account  for 
the  excessive  retreat  of  the  ice  edge  in  summer.  A  higher 
albedo,  produced  by  a  more  "accurate"  albedo  representation, 
would  reduce  the  excessive  melting  and  enhance  freezing.  The 
resulting  increased  ice  thicknesses  would  provide  additional 
resistance  to  earlier  meltoff. 

The  results  of  this  experiment  showed  that  this  was  not 
the  case.  Increasing  the  ice  and  snow  albedo  did  result  in 
nominally  more  ice  area  and  thickness,  especially  near  the 
edges  of  the  pack  ice,  and  less  ice  was  produced  at  reduced 
albedos.  However  the  differences  were  small  even  when  the 
albedos  were  set  well  outside  normal  values. 
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One  explanation  could  be  that  the  upper  ocean  controls  the 
extent  of  the  ice  cover  directly.  With  a  fixed  albedo  of 
0.10,  the  surface  water  absorbs  much  more  heat  than  the  ice. 
Wherever  open  water  and  ice  are  adjacent  during  spring  and 
summer,  the  warmed  water  is  free  to  advect  under  and  around 
the  ice  to  promote  melting  regardless  of  the  albedo  of  the 
ice.  The  reduced  ice  concentration  in  the  MIZ  would  further 
promote  this  process.  However,  this  mechanism  would  also  be 
possible  in  the  earlier  thermodynamic-only  ice  models.  Since 
those  were  the  models  which  displayed  considerable  albedo 
sensitivity,  this  explanation  can  be  rejected. 

A  second  explanation  is  that  the  inclusion  of  the 
interactive  ocean  could  permit  a  feedback  mechanism  which 
largely  compensates  for  any  changes  to  the  frozen-surface 
albedo.  For  example,  the  greater  ice  production  rate  at 
higher  albedos  could  initiate  greater  vertical  mixing  through 
salt  extrusion  and  thereby  bring  more  heat  to  the  surface, 
slowing  the  ice  growth.  Similarly  the  reduced  albedo  would 
initiate  more  melting,  create  a  stronger  halocline  and 
therefore  less  vertical  heat  advection. 

A  third  explanation,  and  probably  the  most  reasonable,  is 
that  the  ocean  heat  flux,  ocean  currents  and  wind  stress  are 
much  more  dominant  in  controlling  the  ice  edge  than  the 
atmospheric  heat  flux  changes  to  the  ice  resulting  from  albedo 
modifications . 
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The  primary  difference  between  this  ice  model,  which  has 
been  shown  to  be  insensitive  to  ice  surface  albedo  changes, 
and  previous  models  which  were  sensitive  to  albedo  changes  was 
the  inclusion  of  a  fully  prognostic  ocean.  Runs  D2  -  D8  were 
used  to  examine  the  impact  of  the  dynamic  features  on  the 
model's  sensitivity  to  frozen-surface  albedo  changes.  The 
greatest  sensitivity  was  observed  in  the  runs  in  which  the 
wind  stress  was  reduced  to  zero  (D6-D8) .  This  suggests  that 
dynamic  forcing,  and,  in  particular  the  wind  stress  might  act 
to  remove,  compensate  or  dominate  most  of  the  effect  of 
changing  the  frozen-surface  albedo.  However  the  model  used 
by  Hibler  (1980)  accounted  for  the  effects  of  wind  stress  in 
simulating  the  annual  cycle  of  Arctic  ice  cover  and  still 
found  that  the  model  was  sensitive  to  small  changes  in  the  ice 
albedo.  The  major  difference  remaining  between  the  Hibler 
model  and  the  model  used  here  was  that  Hibler  used  a  constant 
vertical  ocean  heat  flux.  One  must  therefore  assume  that  it 
is  this  vertical  heat  flux  from  the  ocean,  which  simulations 
from  this  work  indicate  varies  dramatically  in  time  and  space, 
that  reduces  this  model's  sensitivity.  It  follows  then  that 
the  vertical  ocean  heat  flux  is  likely  the  most  dominant 
factor  in  determining  the  thermodynamic  balance  near  the  ice 
edge. 

The  sensitivity  of  runs  D6  and  D7  to  albedo  changes, 
despite  using  a  simulated  variable  ocean  heat  flux,  indicated 
that  another  factor  was  also  important  to  albedo  sensitivity. 
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Ice  thicknesses  and  areal  coverage  were  severely  reduced  in 
run  D6  compared  to  run  D7 ,  and  run  D6  also  displayed  more 
sensitivity  to  the  albedo  change.  The  areas  of  thin  ice 
became  ice  free  much  more  quickly  and  did  not  extend  as  far 
when  the  albedo  was  low.  It  would  therefore  appear  that  the 
model's  albedo  insensitivity  was  also  dependent  on  the  ice 
thickness. 

Further  support  for  this  is  evident  by  comparing  runs  D3 
and  D8 .  In  both  cases  the  ocean  heat  flux  was  reduced  to 
zero;  however  D3  retained  the  dynamic  forcing  from  wind  and 
ocean  currents.  The  lack  of  dynamic  forcing  in  D8  reduced  the 
simulated  ice  thicknesses  in  the  central  pack  and  the 
thickness  gradients  in  the  MIZ.  The  extended,  thin  ice  MIZ 
in  D8  was  much  more  susceptible  to  the  albedo  change  than  the 
strong  thickness  gradient  MIZ  in  D3 . 

There  is  an  apparent  contradiction  in  that  D7  has  thicker 
ice  than  D6  despite  reduced  dynamic  forcing.  However  it 
appears  that  the  water  stress,  which  is  accounted  for  in  Run 
D6  but  not  in  D7 ,  advects  the  ice  away  from  the  thick 
compression  areas  to  regions  of  high  ocean  heat  flux  where  it 
is  melted.  This  reduces  the  total  ice  volume  and  ice 
thicknesses,  especially  in  the  central  pack  and  along  the 
Canadian  boundary  of  the  Beaufort  Gyre.  This  process  may  also 
explain  the  reduced  ice  thicknesses  immediately  adjacent  to 
the  Canadian  Archipelago  and  northern  Greenland  evident  in  the 
ice  strength  experiment  (Chapter  V) . 
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D.   CONCLUSIONS 

The  use  of  only  one  integration  year  vice  several  years  or 
decades  of  integration  obviously  limited  the  degree  of  ice 
cover  modification  possible  between  model  runs  with  different 
surface  albedo  representations.  However  the  objective  here 
was  not  to  determine  what  final  end  condition  of  ice  cover 
could  be  produced  but  whether  or  not  alteration  of  the  albedos 
within  reasonable  limits  could  improve  the  simulation  of 
interannually  varying  ice  cover. 

The  largest  error  in  the  simulated  ice  fields  was  the 
excessive  summer  melt.  A  secondary,  but  still  important 
error,  was  excessive  winter  ice  coverage  in  the  Barents  Sea. 
This  model  appears  quite  insensitive  to  albedo  changes  of  the 
snow/ice  cover.  Surface  albedos  had  to  be  set  well  above 
reasonable  levels  to  produce  even  a  nominal  decrease  in  the 
meltback  in  the  western  Arctic  and  the  same  albedo  changes 
marginally  worsened  the  simulation  of  the  ice  edge  in  the 
Barents  Sea.  Therefore,  changing  the  albedo  representation 
was  not  an  effective  method  of  improving  the  accuracy  of  the 
simulated  ice  cover.  In  fact,  the  insensitivity  of  this  model 
to  albedo  modifications  would  suggest  that  even  the  very 
simple  albedo  representation  initially  used  would  be 
sufficient  and  would  not  degrade  the  model's  accuracy. 

The  primary  reason  for  this  model's  insensitivity  to 
albedo  changes  in  the  frozen  surface  appears  to  be  the 
dominance  of  the  vertical  ocean  heat  flux  in  the  thermodynamic 
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balance  near  the  ice  edge.  Albedo  sensitivity  in  the  model 
is  only  apparent  in  those  few,  small  regions  where  both  the 
vertical  ocean  heat  flux  is  small  and  the  ice  is  thin.  Under 
these  conditions,  the  changes  to  the  thermodynamic  balance 
induced  by  albedo  changes  in  the  ice  cover  are  large  enough 
to  cause  noticeable  changes  in  the  ice  concentration  and 
extent. 

It  should  be  reiterated  here  that  this  study  deals  with 
the  ice  cover  on  an  average  monthly  basis  at  a  fairly  large 
scale.  Albedo  does  have  a  major  impact  on  the  ice  cover  on 
shorter  time  scales  as  evidenced,  for  example,  by  the  flash 
melt  phenomenon.  The  entire  ice  surface  in  a  region  can  melt 
in  a  few  days  when  the  reduced  albedo  of  the  melting  surface 
accelerates  the  melting  process.  However,  on  the  larger  time 
and  space  scales  used  here,  this  type  of  effect  is  largely 
smoothed  out  leaving  the  consistent  longer  term  ocean  heat 
flux  as  the  dominant  controller  of  the  ice  cover. 

The  ocean  currents,  although  not  a  major  factor  in  albedo 
sensitivity,  play  a  significant  role  in  the  evolution  of  ice 
thickness.  They  appear  to  advect  the  ice  away  from  the  thick 
compression  regions  to  regions  of  high  heat  flux  where  the  ice 
is  usually  melted.  This  has  the  effect  of  reducing  the  ice 
thickness  in  the  compression  regions,  particularly  along  the 
northern  limit  of  the  Canadian  Archipelago. 
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VII.   DOMINANT  MECHANISM  EXPERIMENT 

The  albedo  sensitivity  experiment  included  several  runs 
in  which  various  dynamic  and  thermodynamic  processes  were 
removed  from  the  base  linked  model  to  determine  which 
processes  were  primarily  responsible  for  the  insensitivity  of 
the  model  to  albedo  changes.  Although  the  albedo  changes 
continued  to  have  relatively  little  impact  on  the  simulated 
ice  fields,  the  total  ice  field  changed  considerably  from  the 
base  model  output  as  the  different  processes  were  omitted. 

This  experiment  was  conducted  to  determine  which  elements 
of  the  thermodynamic  and  dynamic  forcing  were  dominant  in  the 
linked  ice-ocean  model. 

A.   EXPERIMENTAL  SETUP 

The  same  base  model  as  used  in  the  albedo  sensitivity 
study  was  used  for  this  experiment.  A  single  integration  year 
was  considered  sufficient  due  to  the  high  sensitivity  of  this 
model  to  dynamic  changes.  The  year  1977  was  again  chosen  as 
the  integration  year  because  the  previous  work  had  provided 
a  large  base  of  simulated  1977  ice  field  data  for  comparison. 
The  following  sensitivity  runs  were  examined: 

-  Dl .   Base  model.   10-year  average  ocean.   Starting  from 

Dec  1976  output  fields.    1977  forcing.    Single 
integration  year. 

-  D2 .   Constant  ocean   heat  flux  into   mixed  layer  (2.0  x 

10"^  deg  C/sec,  equivalent  to  approximately  25  W/m^)  . 
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-  D3 .   No  ocean  heat  flux. 

-  D4 .   No  ocean  current  stress  (level  two  (40  m)  ocean 

velocities  =  0.0  m/s) . 

-  D5.   Ocean  heat  flux  =  2.0  x  10*^  and  no  ocean  current 

stress. 

-  D6.   No  wind  stress  (surface  drag  coefficient  =  0.0). 

-  D7 .   No  wind  or  ocean  current  stress. 

-  D8.   No  wind  stress,  no  ocean  current  stress  and  no  ocean 

heat  flux. 

In  cases  D2  through  D8 ,  the  description  for  each  run  indicates 
the  changes  from  the  Base  Model  (run  Dl) . 

B.   RESULTS 

The  winter  (April)  and  summer  (August)  ice  thickness 
contours  from  run  Dl  are  shown  in  Figures  7.1  -  7.2  for 
comparison  with  the  contours  from  runs  D2-D8. 

The  ocean  heat  flux  value  chosen  for  run  D2  was  the 
approximate  median  of  the  simulated  average  annual  flux  cycle 
over  the  Arctic  region.  Applying  this  constant  value  of  25 
W/m^  over  the  entire  grid  heavily  exaggerated  the  basin-wide 
simulated  ice  cycle.  The  winter  ice  cover  expanded  much  too 
far  south  into  the  Norwegian  Sea  and  the  summer  ice  cover  was 
reduced  to  a  small  fraction  of  the  observed  cover  (Figures 
7.3  -  7.4).  The  central  ice  pack  normally  receives  little  or 
no  ocean  heat  flux.  The  relatively  large  heat  flux  applied 
in  this  run  rapidly  reduced  the  pack  ice  thickness  and  areal 
coverage.  By  August,  the  central  pack  was  almost  completely 
melted  off.   The  subsequent  freezing  in  early  winter  produced 
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Figure  7.1   Simulated  ice  thickness  contours  (cm)  for  April, 
1977.   Base  model  (run  Dl)  . 
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Figure  7.2   Simulated  ice  thickness  contours  (cm)  for  August, 
1977.   Base  model  (run  Dl). 
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Figure  7.3  Simulated  ice  thickness  contours  (cm)  for  April, 
1977.  Constant  ocean  heat  flux  =  25  W/m^  (run 
D2)  . 
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Figure  7.4  Simulated  ice  thickness  contours  (cm)  for  August, 
1977.  Constant  ocean  heat  flux  =  25  W/m^  (run 
D2)  .  . 
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a  totally  unrealistic  simulated  winter  ice  cover.  Although 
the  prescribed  heat  flux  in  the  Barents  and  Norwegian  Seas 
during  winter  was  positive,  it  was  still  considerably  less 
than  in  the  base  model  case  (run  Dl)  .  This  resulted  in 
considerably  less  melting  and  extended  ice  coverage  there. 

Run  D3  used  an  ocean  heat  flux  set  to  zero.  This  produced 
even  heavier  ice  during  winter  in  the  North  Atlantic  and 
Greenland  Sea  (Figure  7.5)  because  the  ocean  heat  flux  in  that 
region  was  reduced  even  further  from  the  values  used  in  run 
Dl.  The  ice  also  remained  longer  and  stayed  thicker  around 
the  western  ice  boundary  where  excessive  meltback  occurred  in 
run  Dl  in  summer  (Figure  7.6)  .  However  the  meltback  was  still 
somewhat  severe  compared  to  the  observed  data.  For  the 
purposes  of  these  observations  the  western  ice  boundary 
included  the  Beaufort,  Chukchi  and  East  Siberian  Seas.  The 
central  pack  thickness  and  distribution  remained  relatively 
unchanged  from  the  base  model  as  the  normal  simulated  heat 
flux  there  was  also  zero. 

Run  D4  had  no  stress  from  ocean  surface  currents.  In  this 
case  the  pack  ice  thicknesses  decreased  slightly  but  the 
relative  distribution  of  the  ice  thickness  changed  very 
little.  In  comparison  to  the  base  run,  the  ice  cover  over  the 
Beaufort  Gyre  was  rotated  counter-clockwise  approximately 
fifteen  degrees  (Figures  7.7  -  7.8).  The  MIZ  expanded  outward 
in  all  regions  except  the  EGG  where  the  ice  advection  was 
suppressed.   The  thickness  gradient  marking  the  poleward  side 
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Figure  7.5   Simulated  ice  thickness  contours  (cm)  for  April, 
1977.   Ocean  heat  flux  =  0.0  W/m^  (run  D3). 
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Figure  7.6   Simulated  ice  thickness  contours  (cm)  for  August, 
1977.   Ocean  heat  flux  =  0.0  W/m^  (run  D3). 
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Figure  7.7  Simulated  ice  thickness  contours  (cm)  for  April, 
1977.  Level  two  current  velocities  =  0.0  m/s  (run 
D4)  .  . 
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Figure  7.8   Simulated  ice  thickness  contours  (cm)  for  August, 
1977.   Level  two  current  velocities  =  0.0  m/s  (run 

D4)  . 
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of  the  Transpolar  Drift  Stream  (TPD)  also  became  less  intense. 
Three  supplementary  runs  were  conducted  with  the  same  dynamic 
conditions  as  run  D4  but  using  different  ice  strength  values 
as  in  the  ice  strength  experiment.  As  ice  strength  was 
increased,  which  reduced  the  plasticity  of  the  ice,  the  main 
features  separating  run  D4  from  run  Dl  became  less  noticeable. 

In  run  D5,  a  similar  counter-clockwise  rotation  to  the  ice 
field  as  in  run  D4  was  observed  for  the  Beaufort  Gyre  ice 
cover.  The  reduced  advection  in  the  TPD  and  EGC  was  evidenced 
by  the  reduced  thickness  gradients  in  the  current  regions  (See 
Figure  7.9).  The  combined  reductions  in  ice  pack  thickness 
from  higher  ocean  heat  flux  in  the  central  Arctic  and  no 
dynamic  ocean  forcing  resulted  in  a  complete  ice  meltback  by 
August. 

The  direct  dynamic  forcing  of  the  wind  on  the  frozen 
surface  was  removed  in  run  D6  by  setting  the  ice/air  and 
snow/air  drag  coefficients  to  zero.  Ice  thickness  was  reduced 
by  approximately  50%  within  six  months  and  the  thickness 
distribution  pattern  also  changed  substantially  (Figures  7.10- 
7.11).  The  central  pack  shifted  over  1000  km  towards  the 
Beaufort  Sea,  but  due  to  the  reduced  ice  thickness,  the 
thickness  gradient  there  changed  very  little.  Ice  remained 
in  the  western  boundary  seas  at  least  one  month  longer  than 
the  standard  case  and  the  initial  MIZ  retreat  in  that  region 
was  not  as  severe.  However  the  overall  reduced  thickness 
resulted  in  a  smaller  total  ice  cover  in  the  Arctic  by  August. 
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Figure  7.9  Simulated  ice  thickness  contours  (cm)  for  April, 
1977.  Ocean  heat  flux  and  level  two  velocities 
=  0.0.  (run  D5)  . 
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Figure  7.10  Simulated  ice  thickness  contours  (cm)  for  April, 
1977.   Wind  stress  =  0.0.  (run  D6) . 
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Figure  7.11  Simulated  ice  thickness  contours  (cm)  for  August, 
1977.   Wind  stress  =  0.0.  (run  D6) . 
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The  advection  of  ice  in  the  TPD  and  EGC  was  noticeably 
reduced . 

The  run  D7  simulation  removed  all  direct  dynamic  forcing 
from  both  the  air  and  ocean.  No  horizontal  ice  velocity  was 
permitted.  The  same  large  shift  of  the  central  ice  pack 
towards  the  Beaufort  Sea  was  observed  similar  to  run  D6  but 
the  ice  thicknesses  were  not  as  severely  reduced.  The  MIZ 
retreat  in  the  Beaufort  Sea  started  earlier  but  did  not  move 
as  far  poleward  as  in  run  D6.  Run  D7  showed  a  noticeable  ice 
retreat  poleward  of  Spitsbergen  and  in  the  EGC  which  would 
correspond  to  the  removal  of  the  TPD  and  EGC  (Figures  7.12  - 
7.13)  .  A  comparison  of  the  summer  thickness  contours  for  runs 
D5,  D6  and  D7  shows  that  in  those  simulations  where  the 
central  pack  is  reduced  to  thicknesses  on  the  order  of  150  cm 
or  less,  the  subsequent  evolution  of  the  ice  cover  changes 
dramatically.  There  appears  to  be  a  critical  thickness  below 
which  the  influence  of  the  various  mechanisms  controlling  the 
ice  changes  significantly. 

Finally  run  D8  used  the  same  dynamic  limitations  as  run 
D7 ;  however  the  ocean  heat  flux  was  also  set  to  zero.  The 
same  excessive  ice  cover  over  the  Barents  Sea  as  was  observed 
in  run  D3  was  repeated  and  the  same  general  ice  thickness  and 
distribution  pattern  as  for  run  D7  was  produced.  However  in 
contrast  to  all  the  other  runs,  the  ice  remained  over  the 
western  boundary  seas  for  the  entire  year-long  simulation. 
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Figure  7.12  Simulated  ice  thickness  contours  (cm)  for  April, 
1977.  Wind  and  ocean  current  stress  =  0.0.  (run 
D7)  . 
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Figure  7.13  Simulated  ice  thickness  contours  (cm)  for  August, 
1977.  Wind  and  ocean  current  stress  =  0.0.  (run 
D7)  . 
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Furthermore,  the  ice  edge  did  not  retreat  poleward  of 
Spitsbergen,  a  condition  closer  to  the  observed  fields  than 
the  previous  simulations  (Figures  7.14  -  7.15). 

C.   DISCUSSION 

This  experiment  was  conducted  primarily  to  determine 
whether  ocean  heat  flux  or  water  stress  from  ocean  currents 
was  the  more  important  process  determining  the  limit  and 
thickness  distribution  of  the  ice  cover.  These  are  the  two 
primary  processes  which  become  more  realistic  once  the 
prognostic  ocean  is  included  in  the  linked  model.  The  results 
indicated  that  such  a  simple  determination  was  not  possible. 
As  has  been  noted  in  the  previous  chapters,  there  appear  to 
be  many  factors  working  in  combination  and  often  with  opposing 
effects . 

The  direct  dynamic  influence  of  ocean  surface  currents 
influences  the  simulated  ice  cover  by  distorting  it  in  the 
direction  of  the  surface  currents.  The  Canadian  Basin  ice 
pack  is  rotated  clockwise  by  the  Beaufort  Gyre  and  ice  is 
advected  along  the  streamlines  of  the  TPD  and  the  EGC.  This 
increases  the  thickness  gradient  poleward  of  the  TPD  and 
pushes  the  summer  ice  edge  in  the  eastern  Arctic  tow.irds 
Spitsbergen  and  well  south  along  the  east  coast  of  Greenland. 
The  region  of  thickest  multi-year  ice  is  shifted  towards  the 
Fram  Strait.  The  ocean  currents  converge  the  ice  to  the 
center  of  the  gyre  and  into  a  region   off  the  north  coast  of 
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Figure  7.14  Simulated  ice  thickness  contours  (cm)  for  April, 
1977.  Wind  and  ocean  current  stress  and  ocean 
heat  flux  =  0.0.  (run  D8). 
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Figure  7.15  Simulated  ice  thickness  contours  (cm)  for  August, 
1977.  Wind  and  ocean  current  stress  and  ocean 
heat  flux  =  0.0.  (run  D8). 
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Greenland  and  the  Canadian  Archipelago.  This  accounts  for  at 
least  a  portion  of  the  thicker  ice  evident  in  those  regions. 

The  ice  is  advected  parallel  to  the  coastline  all  around 
the  perimeter  of  the  Canadian  Basin.  This  reduces  the  ice 
thickness  iminediately  poleward  of  the  Canadian  Archipelago 
but  the  ice  does  not  collect  in  the  Beaufort  Sea  since  most 
of  it  is  quickly  melted  by  the  high  ocean  heat  flux  there. 
The  ice  that  is  left  continues  to  be  pushed  along  the 
coastline  out  of  the  Beaufort  and  Chukchi  Seas  and  into  the 
Laptev  Sea. 

The  ice  cover  in  the  Norwegian  and  Barents  Seas  is  nearly 
unaffected  when  the  water  stress  from  ocean  currents  is 
ignored.  However  changes  to  the  ocean  heat  flux  result  in 
major  changes  to  the  simulated  ice  cover.  Of  course  the  heat 
flux  at  any  position  is  a  function  of  earlier  horizontal  and 
vertical  advection.  Therefore  the  currents  do  serve  a  purpose 
in  this  region.  However  it  is  notable  that  the  dynamic  impact 
of  these  currents  at  the  space  and  time  scales  used  here  is 
small . 

The  wind  stress  appears  to  have  a  much  stronger  effect  on 
the  ice  cover  than  the  ocean  current  stress  at  these  spatial 
and  temporal  scales.  The  ice  thickness  and  distribution  is 
drastically  altered  when  wind  stress  is  removed  from  the 
simulations.  The  wind  stress  appears  to  provide  the  majority 
of  the  compressive  forcing  which  produces  the  realistic  ice 
thicknesses  in  the  central  pack  and  convergence  regions.   For 
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the  1977  forcing  year,  the  wind  provided  an  offshore  component 
to  the  ice  edge  position  around  the  entire  western  half  of  the 
Arctic  Basin,  particularly  in  the  East  Siberian  Sea.  In  the 
eastern  Arctic  Basin  the  wind  stress  appeared  to  increase  the 
ice  advection  along  the  TPD  and  the  EGC  and  the  ice  cover  was 
pushed  southward  in  the  Barents  Sea.  It  should  be  noted  that 
the  Walsh  monthly  mean  wind  fields  used  in  this  study  vary 
interannually  in  direction  and  intensity  (Walsh,  private 
communication) .  Therefore  considerable  differences  in  the 
wind  forcing  and  corresponding  changes  to  the  ice  fields  would 
probably  be  observed  had  a  different  forcing  year  been  used 
in  this  experiment.  Nevertheless,  the  relative  strength  of 
this  forcing  mechanism  would  be  expected  to  remain  similar. 

The  ocean  heat  flux,  wind  stress  and  ocean  currents  form 
a  feedback  which  amplifies  the  influence  of  each  of  the 
individual  mechanisms  on  the  melt  cycle.  For  example,  the 
ocean  heat  flux  thins  the  ice  which  then  allows  the  dynamic 
forcing  to  have  more  effect.  Conversely,  the  thicker  ice 
regions  which  develop  over  areas  of  reduced  or  negative  ocean 
heat  flux  become  less  sensitive  to  the  dynamic  forcing.  This 
feedback  is  stronger  when  the  ice  strength  is  reduced  as 
described  in  the  ice  strength  experiment. 

D.   CONCLUSIONS 

The  simulations  conducted  in  this  experiment  support  the 
view   that   the   importance   of   the   various   dynamic   and 
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thermodynamic  mechanisms  controlling  the  monthly  variations 
of  ice  coverage  vary  regionally  across  the  Arctic  Basin.  The 
ocean  heat  flux  appeared  to  be  the  dominant  mechanism 
controlling  the  extent  of  the  monthly  average  ice  cover  while 
the  wind  and  current  stresses  appeared  to  be  the  mechanisms 
determining  the  ice  thickness  and  distribution.  However  the 
dynamic  forcing  from  winds  and  currents  had  quite  different 
effects  in  the  different  regions  of  the  Arctic  Basin.  In  the 
western  basin  which  was  incorporated  in  regions  1  and  2,  the 
wind  and  currents  acted  in  concert  with  the  ocean  heat  flux 
to  clear  the  ice  from  the  boundary  seas.  In  contrast,  these 
same  forces  acted  to  increase  the  ice  cover  in  the  TPD  and  EGC 
(region  3)  while  in  the  Barents  and  Norwegian  Sea  (region  4), 
the  wind  stress  tended  to  oppose  the  effects  of  ocean  heat 
flux  and  current  stress  had  little  impact.  In  general,  the 
forcing  from  ocean  currents  was  usually  weaker  than  from  ocean 
heat  flux  or  wind  stress.  The  exceptions  were  in  the  strong 
current  regions  along  the  coast  in  the  Canadian  Basin  and 
along  the  axis  of  the  TPD  and  EGC.  Inclusion  of  wind  and 
current  stresses  always  tended  to  increase  the  ice  area, 
either  directly  or  indirectly.  In  regions  3  and  4  the  effect 
was  direct  by  off-ice  advection.  In  regions  1  and  2  the 
effect  was  indirect.  The  dynamics  initially  caused  the  ice 
to  converge  which  thinned  the  ice  on  the  boundaries  thereby 
allowing  more  ice  growth  and  subsequent  convergence.  As  the 
central  pack  became  thicker,  it  also  became  stiffer.   The 
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thickness  gradients  began  to  extend  further  from  the  central 
core  and  the  MIZ  also  increased  in  thickness.  The  ocean  flux 
was  unable  to  melt  the  same  extent  of  a  thicker  MIZ;  therefore 
the  ice  area  expanded. 

The  ocean  does  not  appear  to  control  the  ice  coverage 
significantly  through  direct  dynamic  forcing.  However  it  does 
appear  to  have  a  very  important  thermodynamic  role.  The  three 
dimensional  ocean  circulation  positions  the  heat  contained  in 
the  intermediate  depth  ocean  layers  in  the  appropriate 
regions.  This  heat  is  then  available  whenever  the  vertical 
circulation,  driven  by  conditions  at  the  surface,  extends  deep 
enough  to  tap  it.  The  ocean  essentially  preconditions  the 
ocean  below  the  mixed  layer  so  that  the  heat  necessary  to 
control  the  ice  edge  at  any  location  is  available. 

The  dramatic  changes  to  the  ice  cover  evident  in  the 
simulations  where  wind  stress  and  ocean  heat  flux  were 
modified  emphasized  the  sensitivity  of  the  model  to  these 
parameters.  The  importance  of  realistic  ocean  heat  flux  was 
already  suggested  in  previous  chapters  and  is  further 
supported  by  the  results  of  this  experiment.  This  experiment 
also  emphasized  the  importance  of  appropriate  wind  forcing 
fields  if  the  monthly  average  ice  cover  is  to  be  simulated 
accurately  in  the  Arctic. 


188 


VIII.   STATISTICAL  CONSIDERATIONS 

The  time  series  of  "observed"  ice  concentration  data 
covered  the  same  regions  at  the  same  grid  scale  and  over  the 
same  time  period  as  the  simulated  data.  This  permitted 
calculation  of  the  correlations  between  the  various  model 
simulations  and  the  observed  fields.  The  differences  between 
the  correlations  calculated  for  each  new  experiment  were 
compared  to  determine  if  the  changes  incorporated  into  the 
model  improved  the  accuracy  of  the  simulated  ice  concentration 
fields.  A  large  increase  in  the  correlation  coefficient 
indicated  that  the  modifications  made  to  the  model  improved 
the  accuracy  of  the  output.  The  sign  of  the  correlation 
differences  were  usually  consistent  indicating  that  the 
modifications  made  had  produced  a  noticeable  and  consistent 
change.  However  statistical  analysis  of  the  results  was 
required  in  order  to  determine  if  these  changes  were 
statistically  significant. 

Chervin  and  Schneider  (1976)  and  Chervin  (1980) 
investigated  a  similar  problem  as  applied  to  atmospheric 
general  circulation  models  (GCMs) .  They  suggested  that  a 
knowledge  of  the  "noise  climatology"  or  natural  variability 
inherent  in  a  GCM  was  a  prerequisite  for  judging  the 
significance  of  the  results  from  prescribed  change  model 
experiments.  Change  model  experiments  are  similar  to  what  has 
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been  done  in  this  work  in  which  various  elements  of  a  base 

model  are  changed  to  investigate  changes  in  the  simulated  end 

conditions.  Chervin  and  Schneider  (1976)  noted  that  the  basic 

problem  was  one  of 

distinguishing  between  signal  (that  part  of  the  difference 
between  the  results  of  a  prescribed  change  experiment  and 
an  unperturbed  control  case  which  is  attributable  to  the 
prescribed  change)  and  noise  (some  measure  of  inherent 
model  variability)  in  a  prescribed  change  response. 

The  model  variability,  in  the  case  of  the  atmospheric 

GCM's  referred  to  above,  was  a  result  of  the  numerical 

representation  of  random  processes  in  the  atmosphere.   Each 

run  of  such  a  model  allows  for  some  randomness  and  as  a  result 

the  simulated  fields  will  vary  from  run  to  run  with  no  change 

in  boundary  conditions  or  numerical  code  within  the  model. 

The  numerical  model  used  here  does  not  have  a  similar  output 

variability.    The  simulated  ice  fields  produced  will  not 

change  unless  the  boundary  conditions  or  numerical  code  is 

modified.   However  in  this  work,  it  is  not  the  model  output 

fields  which  are  directly  compared  but  the  difference  between 

models  of  the  correlations  between  simulated  and  observed  ice 

concentration  fields.    Therefore  it  is  these  correlation 

differences   which  must   be   examined   to   determine   their 

statistical  significance. 

A.   METHOD 

The  correlations  calculated  in  this  work  use  120  sample 
pairs  (every  month  for  ten  years) .   Each  value  in  each  paired 
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sample  was  the  average  of  several  hundred  grid  values  of  ice 
area  contained  within  each  of  the  four  designated  regions. 
A  common  method  to  establish  the  variability  and  subsequently 
the  confidence  intervals  of  the  correlation  differences, 
required  that  subsets  of  the  paired  samples  be  selected  and 
assumed  independent.  The  time  series  of  120  paired  samples 
were  divided  into  subsets  which  were  deemed  independent  and 
correlations  calculated  for  each.  Means  and  standard 
deviations  for  these  subset  correlations  were  then  determined. 
However  if  the  number  of  subsets  was  too  small  or  if  the 
subsets  were  not  independent  of  each  other,  their  means  were 
likely  to  be  biased.  Significant  persistence  of  ice 
concentration  anomalies  has  been  observed  in  areas  of  the 
eastern  Arctic  to  lags  of  several  months  and  in  exceptional 
cases  to  over  a  year  (Fleming,  1987) .  The  method  presented 
here  strikes  a  compromise  balancing  the  opposing  requirements 
of  independence  versus  a  fairly  large  number  of  subsets.  A 
subset  size  of  one  year  was  selected  because  it  was 
exceptional  to  find  dependence  at  lags  over  a  year.  This 
determined  that  the  maximum  number  of  subsets  was  10. 

The  standard  formulation  used  to  calculate  the  95% 
confidence  interval  was: 


XB  ±  to_975  X  SD  X  (n)^^^ 


where  XB  was  the  average  correlation  over  10  subsets,  t  was 
the  Student  t  distribution  value  at  95%  confidence  level  and 
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9  degrees  of  freedom,  SD  was  the  estimated  standard  deviation 
of  the  10  subset  correlations  and  n  was  the  number  of  years 
(10)  .  This  formula  was  derived  from  normal  theory  which 
required  n  >  20  and  observations  independent  and  identically 
distributed.  Obviously  our  sample  size  of  10  made  this 
technique  suspect,  however  it  was  the  largest  sample  possible 
in  view  of  the  dependence  consideration. 

An  alternative  technique  used  to  investigate  the 
statistical  significance  of  the  difference  between  the  various 
correlations  was  called  "Jackknife."  (For  a  complete 
description  see,  for  example,  Mosteller  and  Tukey  (1977)). 
This  technique  offers  a  way  to  establish  sensible  confidence 
limits  in  complex  situations  in  which  the  observations  are  not 
necessarily  independent  and  identically  distributed.  The 
correlations  were  calculated  with  all  the  data  and  then  after 
dividing  the  data  into  subsets,  correlations  were  recalculated 
by  leaving  out  each  subset,  one  at  a  time.  The  subset  size 
chosen  was  again  one  year.  Pseudo  values  were  determined  by 
weighting  the  correlations  appropriately  to  account  for  the 
sampling  "overlap"  and  the  Student  t  tables  were  then  used  to 
establish  the  confidence  limits.  The  procedure  is  outlined 
below. 

The  correlations  calculated  for  two  different  models  using 
all  10  years  of  data  are  r^(0)  and  r2(0).  The  correlation 
difference,  Y(0),  is  then 

Y(0)  =  ri(0)  -  r2(0) 
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If  r^  (n)  and  r2(n)  are  the  correlations  leaving  out  year  n, 
where  n=l, 2 , . . . , 10,  then 

Y(n)  =  ri(n)  -  r2(n) 

The  n*^  pseudo  value,  Y  (n)  ,  is  calculated  as  follows: 

Y*(n)  =  (10  X  Y(0))  -  (9  X  Y(n)) 

where  10  and  9  are  the  weights  appropriate  to  our  one  year 
subset  size.  The  variances,  s  ,  of  the  ten  pseudo  values  are 
calculated  and  the  approximate  two-sided  95%  confidence  limit 
for  the  difference  of  the  true  correlation  was  given  by: 

Y(0)  ±  to.975  X  (sVlO)'^^ 

where  t^  ^j^  =  2.2  62.  This  value  was  found  in  the  Student  t 
table  using  9  degrees  of  freedom. 

B.   RESULTS 

The  common  and  Jackknife  techniques  were  applied  to  the 
correlations  calculated  in  the  Interactive  Ocean  Experiment 
(Table  5) .  For  the  common  method,  means,  standard  deviations 
and  confidence  intervals  were  calculated  from  each  group  of 
ten  subset  correlations.  These  are  shown  in  Table  8. 
Confidence  intervals  were  calculated  for  both  sets  of 
correlations  (the  B  series  from  the  model  with  the  interactive 
ocean  and  the  C  series  from  the  model  using  the  ten  year 
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TABLE  8 

MEANS,  STD  DEVIATIONS  AND  CONFIDENCE  INTERVALS  OF  THE 

CORRELATIONS  OF  10,  1-YEAR  SUBSETS.   RESULTS  FOR 

EXPERIMENT  B  (FIRST  LINE)  AND  C  (SECOND  LINE) 

IN  ALL  FOUR  REGIONS 


ICE  AREA  TIME  SERIES 


MEANS 

STD  DEV 

95%  CONF.  INT. 

REGION 

CORR  B/CORR  C 

CORR  B/CORR  C 

CORR  B/CORR  C 

1 

.9322 

.0580 

.8907-. 9737 

.9200 

.0504 

.8839-. 9561 

2 

.9323 

.0529 

.8945-. 9701 

.9198 

.0627 

.8750-. 9646 

3 

.9628 

.0164 

.9511-. 9745* 

.9396 

.0271 

.9202-. 9590* 

4 

.9500 

.0190 

.9364-. 9636* 

.9061 

.0466 

.8728-. 9394* 

ANOMALY  TIME  SERIES 


MEANS 

STD  DEV 

95%  CONF.  INT. 

REGION 

CORR  B/CORR  C 

CORR  B/CORR  C 

CORR  B/CORR  C 

1 

.3472 

.4404 

.0322-. 6622 

.2779 

.3864 

.0015-. 5543 

2 

.7198 

.1696 

.5985-. 8411* 

.5670 

.3614 

.3085-. 8255 

3 

.4021 

.2422 

.2289-. 5753 

.3305 

.2583 

.1457-. 5153 

4 

.6266 

.2428 

.4529-. 8002 

.4952 

.2354 

.3268-. 6636 

B  simulation  was  with  fully  interactive  ocean  model 

C  simulation  was  with  10-year  mean  cycle  ocean 

*  difference  between  the  correlations  is  95%  significant 
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average  ocean  condition) .  The  intervals  were  then  compared 
to  determine  if  they  confirmed  each  other  and  how  important 
the  variance  of  the  correlations  was. 

The  differences  between  the  correlations  for  regions  3  and 
4  of  the  ice  area  time  series  were  determined  to  be 
significant  at  the  95%  confidence  level.  The  average  of  the 
C  correlations  fell  outside  the  confidence  limits  of  the  B 
correlations  and  the  average  of  the  B  correlations  fell 
outside  the  confidence  interval  of  the  C  correlations.  This 
provided  some  measure  of  confirmation.  The  correlation 
differences  between  B  and  C  for  the  ice  area  anomaly  time 
series  were  primarily  insignificant  at  the  95%  confidence 
level.  One  exception  was  the  B  model's  region  2.  In  that 
case,  the  C  model  average  correlation  fell  outside  the  B 
model's  correlation  confidence  interval.  However,  the 
opposite  was  not  true  due  to  a  much  higher  correlation 
variance  for  the  C  model. 

In  the  process  of  applying  the  Jackknife  technique,  the 
correlation  differences  (B  subtract  C)  for  all  four  regions, 
for  both  time  series,  and  for  10  pseudo  value  subsets  were 
calculated.  These  are  shown  in  Table  9.  Note  that  they  were 
consistently  positive.  Table  10  shows  the  calculated  pseudo 
values,  their  means,  their  standard  deviations  and  the 
calculated  95%  confidence  intervals  for  the  ice  area  time 
series.   Table  11  shows  the  same  information  for  the  ice  area 
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TABLE  9 

CORRELATION  DIFFERENCES  BETWEEN  B  AND  C  FOR  THE  TOTAL  TIME 
SERIES  AND  FOR  THE  10  PSEUDO  SUBSETS  (ONE  YEAR  REMOVED) 


ICE  AREA 

TIME  SERIES 

PSEUDO 

SUBSET 

REGION  1 

REGION  2 

REGION  3 

REGION  4 

TOTAL 

.0274 

.0464 

.0355 

.0854 

1 

.0295 

.0328 

.0388 

.1079 

2 

.0266 

.0453 

.0290 

.0883 

3 

.0340 

.0494 

.0388 

.0540 

4 

.0221 

.0461 

.0341 

.0941 

5 

.0316 

.0579 

.0353 

.0809 

6 

.0270 

.0469 

.0436 

.0888 

7 

.0306 

.0533 

.0362 

.0928 

8 

.0248 

.0482 

.0306 

.0792 

9 

.0233 

.0306 

.0322 

.0825 

10 

.0233 

.0494 

.0349 

.0814 

ICE  AREA  ANOMALY  TIME  SERIES 


PSEUDO 

SUBSET 

REGION  1 

REGION  2 

REGION  3 

REGION  4 

TOTAL 

.0973 

.1780 

.2336 

.2840 

1 

.1093 

.1038 

.2496 

.3731 

2 

.0922 

.2539 

.1826 

.3110 

3 

.1271 

.1647 

.2911 

.1577 

4 

.0974 

.1978 

.2185 

.3166 

5 

.1407 

.2125 

.2104 

.2498 

6 

.0976 

.2034 

.2978 

.2807 

7 

.1257 

.1956 

.2644 

.3094 

8 

.0819 

.1714 

.1586 

.2004 

9 

.0792 

.1068 

.2240 

.3455 

10 

.0312 

.1797 

.2263 

.2971 
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TABLE  10 

PSEUDO  VALUES  FOR  THE  ICE  AREA  TIME  SERIES 
MEANS,  STD  DEVIATIONS  AND  95%  CONFIDENCE  INTERVALS  OF  THE 

PSEUDO  VALUES 


ICE  AREA  TIME  SERIES 


YEAR 

REMOVED 

REGION  1 

REGION  2 

REGION  3 

REGION  4 

1 

.0085 

.1688 

.0058 

-.1171 

2 

.0346 

.0563 

.0940 

.0593 

3 

-.0320 

.0194 

.0058 

.3680 

4 

.0751 

.0491 

.0481 

.0071 

5 

-.0104 

-.0571 

.0373 

.1259 

6 

.0310 

.0419 

-.0374 

.0548 

7 

-.0014 

-.0157 

.0292 

.0188 

8 

.0508 

.0302 

.0796 

.1412 

9 

.0643 

.1886 

.0652 

.1115 

10 

.0643 

.0194 

.0409 

.1214 

STD  DEV 

.0361 

.0756 

.0388 

.1246 

MEANS 

.0285 

.0501 

.0369 

.0891 

CONFIDENCE 

.0016  TO 

-.0077  TO 

.0077  TO 

-.0037  TO 

INTERVAL 

.0532  * 

.1005 

.0633  * 

.1745 

*  indicates  that  the  correlation  difference  is  95%  significant 
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TABLE  11 

PSEUDO  VALUES  FOR  THE  ICE  AREA  ANOMALY  TIME  SERIES 
MEANS,  STD  DEVIATIONS  AND  95%  CONFIDENCE  INTERVALS  OF  THE 

PSEUDO  VALUES 


ICE  AREA  ANOMALY  TIME  SERIES 


YEAR 

REMOVED 

REGION  1 

REGION  2 

REGION  3 

REGION  4 

1 

-.0107 

.8458 

.0896 

-.5179 

2 

.1432 

-.5051 

.6926 

.0410 

3 

-.1709 

.2977 

-.2839 

1.4207 

4 

.0964 

-.0002 

.3695 

-.0094 

5 

-.2933 

-.1325 

.4424 

.5918 

6 

.0946 

-.0506 

-.3442 

.3137 

7 

-.1583 

.0196 

-.0436 

.0554 

8 

.2359 

.2374 

.9086 

1.0364 

9 

.2602 

.8188 

.3200 

-.2695 

10 

.6922 

.1627 

.2993 

.1661 

STD  DEV 

.2794 

.4145 

.3997 

.5883 

MEANS 

.0889 

.1694 

.2450 

.2828 

CONFIDENCE 

-.1026  TO 

-.1185  TO 

-.0523  TO 

-.1368  TO 

INTERVAL 

.2972 

.4745 

.5195 

.7048 
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anomaly  time  series.  The  correlations  are  considered 
significantly  different  using  the  Jackknife  method  if  the 
confidence  interval  does  not  include  zero.  This  is  the  case 
for  region  1  and  3  of  the  ice  area  time  series;  however,  no 
similar  cases  exist  for  the  anomaly  time  series. 

C.   DISCUSSION 

The  consistently  positive  differences  between  the 
correlations  from  the  B  and  C  models  would  suggest  that  the 
B  model  was  a  significant  improvement  over  the  C  model. 
However  neither  of  the  two  statistical  approaches  used  here 
was  able  to  confirm  this  at  the  95%  confidence  level  for  more 
than  three  of  the  eight  correlation  difference  cases.  The 
main  reason  for  this  was  the  high  degree  of  correlation 
variability. 

Correlation  variability  was  most  noticeable  in  the  ice 
area  anomaly  time  series  because  the  dominating  influence  of 
the  annual  cycle  was  removed,  thereby  reducing  the  correlation 
coefficient  values  and  allowing  more  correlation  variability. 
This  was  also  apparent  by  comparing  the  actual  time  series 
plots  (Figures  4.19  -  4.22  with  4.31  -  4.38)  .  The  variability 
indicates  that  the  model  did  very  well  for  some  years  in  some 
regions  yet  very  poorly  in  others.  The  simulated  ice 
condition  at  any  time  was  a  direct  function  of  the  forcing  and 
initial  conditions.  The  model  had  no  inherent  variability. 
Therefore  the  changing  performance  of  the  model  had  to  be  a 
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result  of  some  other  cause.  One  explanation  could  be  that 
some  mechanism  which  was  not  properly  represented  in  the  model 
was  controlling  the  ice  and  varying  dramatically  in  importance 
from  year  to  year,  and  region  to  region.  The  model  would 
perform  fine  when  this  mechanism  had  little  influence  but 
would  be  unable  to  account  for  those  periods  or  regions  when 
the  mechanism  was  important.  Alternatively,  the  prescribed 
forcing  may  be  in  error.  In  particular  the  monthly  averaged 
atmospheric  forcing  may  have  too  long  a  temporal  resolution 
to  adequately  represent  strong,  short  term  episodic  changes. 
The  ice  dynamics  experiment  indicated  that  these  could  be 
important  because  wind  stress  was  able  to  alter  the  ice  field 
dramatically  in  short  periods  of  time  and  the  subsequent 
evolution  of  the  ice  field  was  also  modified.  Of  course  some 
combination  of  model  limitation  and  forcing  error  was  also 
possible. 

D.   CONCLUSIONS 

This  work  was  done  to  determine  if  the  correlation 
differences  between  the  different  model  simulations  and 
observed  fields  of  ice  concentration  could  be  shown  to  be 
statistically  significant.  In  general,  it  appears  that  the 
correlations  determined  between  modelled  and  observed  fields 
of  ice  area  are  apparently  not  yet  stable  enough  to  do  so. 
The  improvement  in  the  correlations  from  model  C  to  model  B 
was  consistent  in  all  regions,  for  all  subsets  and  for  both 
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the  ice  area  and  anomalous  ice  area  time  series.  This  would 
certainly  suggest  that  the  improvement  was  real.  However  the 
correlations  were  so  variable,  particularly  in  the  ice  area 
anomaly  case,  that  supporting  statistics  could  not  be 
provided. 
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IX.   SUMMARY 

This  work  has  involved  the  development  of  a  linked  ice- 
ocean  numerical  model  capable  of  simulating  the  annual  cycle 
and  interannual  variations  of  ice  cover  in  the  Arctic.  The 
accuracy  of  the  simulated  seasonal  ice  area  and  ice  edge 
position  appears  to  be  an  improvement  over  any  other  large 
scale  models  operating  at  present,  yet  the  computational 
requirements  have  been  maintained  at  very  reasonable  levels. 

A.   MAIN  CONCLUSIONS 

""  '  The  objectives  for  this  thesis  have  been  met.  It  has  been 
shown  that  the  inclusion  of  an  interactive,  prognostic  ocean 
component  in  the  ice  model  provides  substantial  improvement 
to  simulations  of  both  the  annual  cycle  and  interannual 
variations  of  the  ice  cover.  Closer  examination  of  the 
effects  of  the  various  ice  control  mechanisms  indicated  that 
the  main  reason  for  improved  ice  cover  simulation  with  the 
prognostic  ocean  was  a  more  realistic  representation  of  the 
ocean  heat  flux.  The  ocean  heat  flux  varied  considerably 
between  regions  and  between  seasons.  The  surface  currents 
were  also  highly  variable.  However,  in  general,  the  forcing 
from  ocean  currents  has  a  second  order  effect  on  the  ice 
cover.  Exceptions  to  this  were  evident  in  the  strongest 
current  regions  (East  Greenland  Current,  Transpolar  Drift 
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stream  and  boundaries  of  the  Beaufort  Gyre)  where  the  dynamic 
impact  of  the  currents  had  a  more  pronounced  effect  on  the 
ice. 

The  linked  model  was  sensitive  to  changes  in  the  ice 
rheology.  Based  on  comparison  with  highly  averaged  observed 
thickness  data  (Bourke  and  Garret,  1987)  the  optimum  value 
chosen  for  the  strength  parameter  in  the  "bulk  viscous" 
rheology  was  P*  =  h  x  10^.  This  value  reduced  the  rigidity  of 
the  ice  pack,  allowing  greater  compression  and  therefore  more 
realistic  ice  thickness  and  distribution.  However,  the 
reduction  in  the  ice  strength  degraded  the  correlation  between 
simulated  and  observed  total  ice  cover.  The  ice  retreat  in 
the  western  Arctic,  which  was  already  somewhat  severe  at  the 
higher  ice  strength,  receded  further  when  P  was  reduced. 
Since  there  was  no  compensating  improvement  in  the  winter  ice 
extent,  the  correlation  coefficients  decreased. 

The  increased  melting  in  the  western  Arctic  at  the 
reduced  ice  strength  appeared  to  be  a  result  of  a  feedback 
between  ice  formation  rates  and  ocean  heat  flux.  The  reduced 
ice  strength  permitted  a  stronger  response  of  the  ice  to 
dynamic  forcing.  This  forcing  was  generally  off-shore  in  the 
western  Arctic.  The  ice  was  forced  to  converge  in  the  central 
pack  and  along  the  north  shore  of  the  Canadian  archipelago, 
opening  up  large  areas  of  open  water  and  thin  ice  over  the 
shelves  in  the  western  Arctic  seas.  This  permitted  faster 
cooling  of  the  surface  water,  greater  ice  production,  more 
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salt  extrusion  and  greater  vertical  mixing.  The  mixing  was 
strong  enough  to  break  through  the  halocline  and  increase  the 
ocean  heat  flux  to  the  surface.  This  resulted  in  more 
extensive  melting. 

The  fully  linked  model  with  the  prognostic  ocean  was  quite 
insensitive  to  albedo  changes  in  the  ice  cover.  This  result 
was  a  clear  contrast  to  previous  model  studies  which 
demonstrated  strong  sensitivity  to  this  parameter.  The  main 
mechanism  included  in  this  model  but  not  in  the  models  which 
displayed  albedo  sensitivity  was  the  interannually  variable 
ocean  heat  flux.  This  term  appears  to  dominate  the 
thermodynamic  balance  near  the  ice  edge.  The  ocean  heat  flux 
thereby  controls  the  concentration  of  the  ice  in  this  region 
and  the  overall  extent  of  the  ice  pack.  The  thermodynamic 
impact  of  changing  surface  albedo  is  relegated  to  minor 
importance . 

Results  from  the  dynamic  mechanism  experiment  supported 
the  idea  that  the  ocean  heat  flux  dominated  the  albedo  effect. 
However  these  results  also  indicated  that  there  were 
conditions,  particularly  when  the  ice  was  thin,  when  the 
albedo  did  become  important.  Those  cases  which  showed  some 
albedo  sensitivity  had  a  combination  of  a  small  ocean  heat 
flux  and  thin  ice.  Under  these  conditions,  the  change  in  the 
thermodynamic  balance  due  to  a  change  in  the  surface  albedo 
could  become  relatively  large  and  the  model's  albedo 
sensitivity  was  increased.   These  conditions  were  possible 
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over  the  shallow  shelves,  in  newly  frozen  leads  within  the 
central  pack  and  along  some  areas  of  the  MIZ.  Realistic 
representation  of  the  frozen-surface  albedo  was  therefore 
considered  worthwhile  and  the  improved  albedo  scheme  of  Ross 
and  Walsh  (1987)  was  incorporated  into  the  model. 

Examination  of  the  relative  importance  of  the  various 
dynamic  and  thermodynamic  ice  cover  forcing  mechanisms 
indicated  that  each  part  of  the  fully  linked  model  had  a 
unique  but  interactively  important  role  in  the  evolution  of 
the  ice  cover.  The  ocean  heat  flux  appeared  to  be  the  overall 
dominant  factor  controlling  the  position  of  the  ice  edge  and 
the  extent  of  the  ice  cover.  Within  the  polar  pack,  it  was 
the  dynamic  forcing,  and  in  particular,  the  wind  forcing  which 
controlled  the  ice  thickness  and  thickness  distribution. 

The  dynamic  portion  of  the  ice  model  can  be  viewed  as 
working  in  concert  with  the  ocean  model  to  produce  the  desired 
ice  cover.  The  ocean  model  does  not  provide  much  of  a  dynamic 
effect  from  surface  currents  but  it  does  provide  an  important 
thermodynamic  control.  The  ocean  circulation  below  the  mixed 
layer  acts  to  position  heat  below  those  regions  where 
observations  indicate  that  the  ice  is  usually  thin  or  melted 
away.  The  dynamic  ice  model  in  turn  tends  to  compress  the  ice 
in  those  areas  where  it  should  regionally  be  thickest  and  thin 
the  ice  in  those  areas  where  the  heat  flux  will  be  used  to 
melt  back  the  ice  cover.    The  linkage  between  these  two 
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processes  is  controlled  by  conditions  at  the  surface  and  the 
response  of  the  mixed  layer. 

Vertical  mixing  needs  to  be  induced  to  tap  the 
intermediate  layer  heat  source  while  stratification  must  be 
induced  to  shut  the  heat  off.  However,  the  large  horizontal 
diffusion  coefficient  required  for  numerical  stability 
probably  causes  excessive  lateral  mixing  of  heat  within  the 
intermediate  layer.  The  increased  heat  available  in  some 
regions  may  be  sufficient  to  maintain  excessive  heat  flux  into 
the  mixed  layer,  despite  stratification  from  fresh  meltwater. 

It  must  be  noted  that  this  assessment  of  the  relative 
importance  of  the  various  mechanisms  is  limited  in  application 
to  the  space  and  time  scales  associated  with  this  model. 
However  Ikeda  et  al.  (1988)  have  presented  similar  conclusions 
based  on  a  model  of  much  smaller  scale  in  the  Labrador  Sea. 
Additionally,  Walsh  et  al.  (1985)  and  Hibler  and  Bryan  (1987) 
also  present  conclusions  consistent  with  those  here  regarding 
the  dominance  of  thermodynamic  processes  in  the  determination 
of  ice  cover  extent. 

Correlations  between  simulated  ice  area  and  observed  ice 
area  were  calculated  for  each  major  change  in  the  model 
configuration.  Two  statistical  approaches  were  used  to 
determine  if  the  difference  between  correlations  obtained  from 
each  different  model  configuration  were  large  enough  to  be 
statistically  significant.  Limited  success  was  achieved  in 
that  the  difference  between  the  two  model  runs  from  the 
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interactive  ocean  experiment,  in  region  3,  were  95% 
significant  for  both  statistical  approaches.  However, 
statistical  significance  could  not  be  determined  at  the  95% 
confidence  level  for  the  majority  of  regions  because  the 
variability  of  the  correlations  was  too  high.  The  large 
correlation  variability  indicated  that  the  model  was  doing 
very  well  in  some  time  periods  but  poorly  in  others.  Two 
explanations  were  proposed.  First,  that  some  process 
controlling  the  ice  area,  for  example  the  mixed  layer,  could 
be  inadequately  represented  by  the  model.  In  those  conditions 
where  the  process  was  relatively  unimportant,  the  model  would 
do  fine.  However,  if  conditions  changed  wherein  the  process 
became  important,  the  model  simulations  would  be  degraded. 
A  second  possible  explanation  is  that  the  atmospheric  forcing 
is  not  being  handled  correctly.  In  particular,  the  monthly 
averaged  wind  stress  is  probably  not  accounting  for  the 
effects  of  short  term,  intense  weather  features.  The  dynamic 
mechanism  experiment  pointed  out  the  large  impact  that  the 
wind  stress  has  on  the  ice  field;  therefore,  this 
underestimation  of  the  wind  stress  could  result  in  periods  of 
poor  model  performance. 

This  work  has  shown  that  the  relative  importance  of  the 
various  processes  controlling  the  ice  cover  in  the  Arctic  is 
regionally  dependent.  For  example,  in  the  western  Arctic  the 
dynamic  effects  assist  the  ocean  heat  flux  in  clearing  the 
western  boundary  seas  of  ice.   In  contrast,  the  dynamics  in 
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the  eastern  Arctic  tend  to  oppose  the  clearing  of  ice  from 
that  region.  That  is,  the  dynamic  forcing  pushes  the  ice  edge 
southward  while  the  ocean  heat  flux  tries  to  melt  the  MIZ 
northward.  Various  other  combinations  of  effects  are  also 
evident  in  smaller  sub-regions  within  the  Arctic  Basin.  In 
order  to  allow  for  these  variations,  an  ice  model  applied  to 
the  Arctic  must  be  able  to  represent  each  mechanism  reasonably 
well  and  be  sensitive  to  the  regional  differences  that  change 
the  importance  of  those  mechanisms.  Inclusion  of  a  fully 
interactive  prognostic  ocean  model  in  a  linked  ice-ocean  model 
is  necessary  to  do  this. 

B.   MODEL  LIMITATIONS 

The  discussion  above  already  pointed  out  a  potential 
problem  if  the  monthly  averaged  wind  fields  did  not  provide 
sufficient  temporal  resolution  of  the  wind  to  account  for 
strong,  episodic  wind  anomalies.  The  limitations  of  the 
forcing  data  must  also  be  considered.  This  must  include  not 
only  the  atmospheric  forcing,  which  is  known  to  have  some 
weaknesses,  but  also  the  prescribed  inflow  at  the  boundaries. 

A  thermodynamic  role  for  the  ocean  has  been  proposed 
wherein  the  ocean  acts  to  precondition  the  waters  below  the 
surface  layer  in  certain  regions  with  heat.  However,  despite 
the  recent  reassessment  of  the  Intermediate  Water  circulation 
by  Aagaard  (1988)  which  would  appear  to  support  the  large 
scale   circulation   simulated   by   the   model   used   here, 
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considerable   ambiguity   on   the   smaller   scale   remains. 
Furthermore,    although   the   simulations   indicate   that 
considerable   interannual   variability   exists   for   the 
geostrophic  surface  currents,  the  accuracy  of  this  variability 
has  not  been  determined. 

The  representation  of  the  mixed  layer,  which  is  the 
critical  link  between  the  ice  and  ocean  portions  of  the  model, 
was  very  simplified.  It  was  noted  above  that  the  processes 
which  occur  above  and  below  the  mixed  layer  generally  act  to 
position  the  heat  under  areas  of  reduced  ice  thickness. 
However  it  is  the  action  of  the  mixed  layer,  largely  driven 
from  the  surface,  which  actually  connects  the  ocean  with  the 
ice.  The  simulation  of  any  such  connection  is  currently 
limited  by  the  vertical  resolution  of  the  upper  ocean  layers 
and  the  simple  representation  of  the  mixing  which  occurs 
there. 

The  drag  coefficients  used  here  are  the  same  for  air/ice 
and  air/water  and  constant  for  the  ice/water.  Recent  studies 
in  the  Arctic  have  shown  that  these  drag  coefficients  vary 
dramatically  with  the  surface  type,  ice  concentration,  ridge 
concentration  and  several  other  factors  (e.g..  Guest  and 
Davidson,  1987) .  The  drag  coefficients  determine  the  stress 
and  are  therefore  important  to  the  evolution  of  the  ice  field. 
Their  simplification  reduces  the  dynamic  complexity  of  the 
forcing. 
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C.   PROPOSED  FUTURE  WORK 

The  surface  mixed  layer  is  the  link  between  the  ice  and 
ocean  sub-models,  and  consequently  has  an  effect  on  the 
majority  of  the  ice  forcing  parameters.  The  relative 
importance  of  many  of  these  parameters  has  been  examined  in 
this  work.  The  most  active  times  of  large-scale  ice  edge 
advance  and  retreat  occur  at  the  start  of  the  freezing  season, 
and  the  middle  of  the  melt  season,  when  the  pack  begins  to 
break  apart.  It  is  also  the  same  times  when  the  mixed  layer 
can  undergo  significant  changes  due  to  convective  and  dynamic 
mixing.  Several  authors  have  noted  the  importance  of 
incorporating  a  realistic  mixed  layer  model  in  model 
simulations  of  sea  surface  temperature  (SST) .  Fleming  (1987) 
and  Garcia  (1988)  have  shown  the  strong  correlation  between 
SST  and  ice  concentration  in  the  Arctic.  It  would  therefore 
be  worthwhile  to  examine  the  importance  of  including  a  mixed 
layer  model  into  the  linked  ice-ocean  model  used  here. 

The  role  of  atmospheric  forcing  in  the  interannual 
variability  of  Arctic  sea  ice  was  examined  in  Hibler  and  Walsh 
(1982)  .  That  paper  concluded  that  simulations  of  interannual 
fluctuations  in  the  ice  cover  were  improved  with  the  inclusion 
of  interannually-varying,  atmospherically-forced  dynamics. 
A  similar  conclusion  was  drawn  by  Walsh  et  al.,  (1985)  using 
a  more  elaborate  model;  however  neither  of  these  efforts 
included  interannual  variations  of  the  ocean.  A  re-analysis 
of  the  role  of  interannually  variable  atmospheric  forcing 
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would  be  valuable  using  this  model,  as  it  includes  an 
interannually  varying  ocean.  Further,  mean  monthly  vice  daily 
atmospheric  forcing  could  be  examined. 

The  drag  coefficients  currently  used  are  quite  simple,  as 
noted  above.  A  sensitivity  study  examining  the  importance  of 
the  drag  coefficients  to  simulations  of  the  ice  cover  would 
be  valuable  and,  if  warranted,  a  more  elaborate  representation 
which  accounts  for  the  coefficient  dependencies  should  be 
developed  and  incorporated  into  the  model.  It  is  anticipated 
that  this  would  generate  more  leads  in  the  ice,  encouraging 
increased  ice  growth  and  a  thicker  ice  cover.  Greater 
compression  might  also  occur  in  the  MIZ  and  rough  surface 
compression  regions  where  the  drag  coefficients  are  larger. 

The  work  proposed  above  would  progress  our  knowledge  of 
ice  cover  simulation  and  Arctic  Ocean  circulation  on  the 
seasonal  time  scale.  The  next  logical  step  would  be 
incorporation  of  this  model  into  an  atmospheric  GCM.  This 
would  be  feasible  using  the  current  code  because  it  has  a 
relatively  small  computation  requirement.  The  computation 
requirements  of  such  a  coupled  GCM  are  expected  to  be  well 
within  the  capabilities  of  the  next  generation  of 
supercomputers. 

Finally,  in  order  to  improve  the  simulation  and  prediction 
of  Arctic  ice  on  a  daily  basis  (a  capability  desired  by  most 
of  the  world's  Navies),  the  effects  of  mesoscale  ocean 
features  and  short  duration  atmospheric  forcing  must  be 
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examined.  To  this  end,  the  next  major  improvement  of  this 
model  should  include  an  increase  in  the  horizontal  resolution 
to  at  least  1/8  or  1/10  degree,  bi-harmonic  diffusion,  an 
increase  in  the  number  of  vertical  levels  to  approximately  40 
and  daily  if  not  six  hourly  atmospheric  forcing. 

D.   FINAL  GEM 

Walsh  et  al.,  (1985)  stated: 

. . .model-derived  trends  (of  Arctic  ice  cover)  may  be 
misleading  in  the  absence  of  a  realistic  treatment  of  ice 
dynamics. 

To  this  I  would  add: 

and  the  direct  and  interactive  effects  of  three- 
dimensional  ocean  circulation. 
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